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ABSTRACT 


The  computer  program  VBREAK  is  developed  to  predict  the  time-dependent, 
two-dimensional  velocity  field  under  normally  incident  breaking  waves  on  beaches 
and  coastal  structures.  To  reduce  computation  time  considerably,  use  is  made  of 
the  depth-integrated  continuity  and  horizontal  momentum  equations.  The  momen¬ 
tum  equation  includes  the  momentum  flux  correction  due  to  the  vertical  variation 
of  the  horizontal  velocity.  The  bottom  shear  stress  is  expressed  in  terms  of  the 
near-bottom  horizontal  velocity  immediately  outside  the  thin  wave  boundary  layer. 
The  third  equation  for  the  momentum  flux  correction  is  derived  from  the  depth- 
integrated  wave  energy  equation.  In  order  to  express  these  three  one- dimensional, 
time-dependent  equations  in  terms  of  the  three  unknown  variables  of  the  water 
depth,  depth-averaged  horizontal  velocity,  and  near-bottom  horizontal  velocity,  the 
normalized  vertical  profile  of  the  horizontal  velocity  is  assumed  to  be  cubic  on  the 
analogy  between  turbulent  bores  and  hydraulic  jumps.  Furthermore,  the  turbulent 
shear  stress  is  assumed  to  be  expressed  using  the  turbulent  eddy  viscosity  whose 
mixing  length  is  projiortional  to  the  water  depth. 

The  three  governing  equations  are  solved  using  the  MacCormack  finite  differ¬ 
ence  method  for  its  simplicity  and  success  in  the  computation  of  hydraulic  jumps. 
The  seaward  and  landward  boundary  algorithms  are  extensions  of  those  used  in  the 
previous  one-dimensional  models  such  as  IBREAK.  The  developed  model  VBREAK 
can  function  as  a  one-dimensional  model  if  the  momentum  flux  correction  term  is 
assigned  to  be  zero  at  the  seaward  boundary.  To  assess  the  applicability  of  the 
MacCormack  method  to  the  shallow-water  continuity  and  momentum  equations. 


VBREAK  is  reduced  to  a  one-dimensional  model  and  compared  with  the  previ¬ 
ously  developed  model  IBREAK.  The  differences  in  the  computed  depth  averaged 
velocity  and  free  surface  are  minimal.  The  MacCormack  method  is  hence  consid¬ 
ered  to  be  efficient  and  accurate  in  the  solution  of  the  shallow- water  continuity  and 
momentum  equations. 

The  quasi  two-dimensional  model  VBREAK  is  compared  with  two  laboratory 
data  sets.  The  first  is  from  the  experiment  conducted  by  Stive  and  outlined  in  Stive 
(1980)  and  Stive  and  Wind  (1982).  Finally,  the  model  is  compared  with  the  data 
collected  and  analyzed  by  Cox  et  al.  (1995).  VBREAK  predicts  the  free  surface  and 
depth  averaged  velocity  well.  The  addition  of  the  momentum  flux  correction  in  the 
governing  equations  has  little  effect  on  these  quantities  as  anticipated  by  Kobayashi 
and  Wurjanto  (1992).  However,  the  vertical  variation  of  the  horizontal  velocity 
is  sensitive  to  this  momentum  flux  correction  because  there  would  be  no  vertical 
variation  without  this  correction  term.  The  phase  speed  is  predicted  reasonably 
well.  The  agreement  between  the  measured  and  computed  vertical  variation  of  the 
horizontal  velocity  is  effected  by  the  mismatch  of  the  measured  and  computed  phase 
speed.  Although  the  energy  dissipation  due  to  wave  breaking  is  modeled  explicitly 
in  VBREAK,  energy  dissipation  in  the  model  is  primarily  numerical  for  breaking 
waves  on  gentle  slopes.  This  is  likely  due  to  the  fact  that  the  energy  dissipation 
and  associated  landward  mass  flux  in  the  surface  roller  are  not  accounted  for  in 
VBREAK.  Nevertheless,  VBREAK  predicts  the  vertical  variation  of  the  horizontal 
velocity  measured  below  the  trough  reasonably  well. 

Originally,  this  work  was  submitted  by  Bradley  D.  Johnson  in  the  Spring  of 
1996  as  a  thesis  (Johnson  1996)  at  the  University  of  Delaware  in  partial  fulfillment 
of  the  requirements  for  the  degree  of  Master  of  Civil  Engineering.  A  summary 
of  this  report  is  being  presented  at  the  25*^  International  Conference  of  Coastal 
Engineering( Johnson  et  al.  1996). 


Chapter  1 


INTRODUCTION 


Available  time-dependent,  one-dimensional  and  other  numerical  models  for 
breaking  and  nonbreaking  waves  on  inclined  structures  and  beaches  were  reviewed 
by  Kobayashi  and  Poff  (1994).  The  one-dimensional  shallow- water  models  (e.g., 
Kobayashi  and  Wurjanto  1989;  Kobayashi  and  Poff  1994)  are  relatively  simple 
and  robust.  Generally,  these  models  predict  the  free  surface  elevation  fairly  ac¬ 
curately,  within  about  20%  errors.  Raubenheimer  et  al.  (1995)  showed  that  the 
one-dimensional,  shallow-water  model  was  in  good  agreement  with  the  variations  of 
wave  spectra  and  shapes  (e.g.,  wave  skewness)  measured  across  the  inner  surf  and 
swash  zones  on  a  gently  sloping  natural  beach. 

The  one- dimensional  models  predict  only  the  depth-averaged  horizontal  ve¬ 
locity.  The  vertical  velocity  may  be  estimated  using  the  two-dimensional  continuity 
equation  together  with  the  computed  depth-averaged  velocity,  while  the  bottom 
shear  stress  may  be  expressed  by  a  quadratic  friction  equation  based  on  the  depth- 
averaged  velocity.  The  comparisons  with  the  experiment  for  regular  waves  spilling 
on  a  rough,  impermeable  1:35  slope  conducted  by  Cox  et  al.  (1995)  indicated  that 
the  horizontal  velocity  measured  below  the  wave  trough  level  was  represented  by  the 
computed  depth-averaged  velocity  reasonably  well.  The  computed  vertical  velocity 
represented  the  measured  vertical  velocity  at  least  qualitatively  except  under  the 
wave  crest.  The  temporal  variation  of  the  bottom  shear  stress  was  predicted  poorly 
because  errors  in  the  computed  horizontal  velocity  were  magnified  in  the  computed 


bottom  shear  stress  and  because  the  bottom  friction  factor  was  not  really  constant. 
These  limited  comparisons  suggest  that  a  vertically  two-dimensional  model  will  be 
required  to  predict  the  detailed  vertical  variations  of  the  fluid  velocities  and  shear 
stress  which  are  essential  for  predicting  cross-shore  sediment  transport  on  beaches 
and  hydrodynamic  forces  acting  on  armor  units  on  coastal  structures  (e.g.,  Tprum 
1994). 

A  simplified  two-dimensional  model  is  developed  in  this  report.  To  reduce 
computational  efforts  considerably,  the  normalized  vertical  variation  of  the  horizon¬ 
tal  velocity  outside  the  wave  boundary  layer  is  assumed  to  be  cubic.  The  vertically 
two-dimensional  problem  is  then  reduced  to  a  depth-integrated  one-dimensional 
problem  in  which  the  three  time-dependent,  one- dimensional  differential  equations 
for  the  water  depth  h,  depth- averaged  horizontal  velocity  U  and  near-bottom  hor¬ 
izontal  velocity  Ub  need  to  be  solved  numerically.  The  simplified  two-dimensional 
model  called  VBREAK  is  computationally  as  efficient  as  the  previous  one-dimensional 
models  such  as  RBREAK2  (Kobayashi  and  Poff  1994).  As  a  result,  VBREAK  can 
be  applied  easily  and  routinely  using  workstations. 

1.1  Outline  of  Report 

The  approximate  governing  equations  adopted  for  VBREAK  are  derived  in 
Chapter  2.  First,  approximate  two-dimensional  equations  for  shallow-water  waves 
on  relatively  gentle  slopes  are  derived  from  the  continuity  and  Reynolds  equations. 
The  approximate  tw'o-dimensional  equations  are  then  integrated  vertically  to  ob¬ 
tain  the  depth-integrated  continuity  and  horizontal  momentum  equations.  This 
momentum  equation  includes  the  unknown  momentum  flux  correction  m  due  to  the 
vertical  variation  of  the  horizontal  velocity  u.  An  equation  for  the  momentum  flux 
correction  m  is  derived  from  the  depth-integrated  wave  energy  equation  (Kobayashi 
and  Wurjanto  1992).  The  bottom  shear  stress  and  wave  energy  dissipation  inside 
the  thin  wave  boundary  layer  are  expressed  in  terms  of  the  near-bottom  horizontal 
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velocity  Uh  and  the  wave  friction  factor  (Jonsson  1966;  Cox  et  al.  1995).  The  ver¬ 
tical  variation  of  the  normalized  horizontal  velocity  u  outside  the  wave  boundary 
layer  is  assumed  to  be  cubic  on  the  analogy  between  turbulent  bores  and  hydraulic 
jumps  (Madsen  and  Svendsen  1983;  Svendsen  and  Madsen  1984).  The  momentum 
flux  correction  m  and  the  wave  energy  dissipation  rate  outside  the  wave  boundary 
layer  due  to  wave  breaking  are  then  expressed  in  terms  oi  h,  U  and  ui.  The  three 
depth-integrated  continuity,  horizontal  momentum,  and  momentum  flux  correction 
equations  may  thus  be  solved  numerically  to  obtain  the  temporal  and  cross-shore 
variations  of  h,  U  and  uj. 

The  numerical  procedures  adopted  to  solve  the  three  governing  equations 
with  appropriate  initial  and  boundary  conditions  are  explained  in  detail  in  Chapter 
3  of  this  report.  The  MacCormack  finite  difference  method  (MacCormack  1969) 
is  selected  because  of  its  simplicity  and  success  in  the  computation  of  unsteady 
open  channel  flows  with  hydraulic  jumps  (Chaudhry  1993).  The  computation  is 
initiated  at  the  time  t  =  0  when  the  specified  incident  wave  train  arrives  at  the 
seaward  boundary  and  no  wave  action  exists  in  the  computation  domain.  The  inter¬ 
val  At  of  each  time  step  for  the  time-marching  computation  is  calculated  using  an 
approximate  stability  criterion  of  the  adopted  explicit  finite  difference  method.  Ap¬ 
proximate  seaward  boundary  conditions  are  used  to  compute  the  boundary  values 
of  h,  U  and  Ub  as  well  as  the  reflected  wave  train  using  the  method  of  character¬ 
istics  (Kobayashi  et  al.  1987,  1989).  The  landward  boundary  algorithm  used  in 
RBREAK2  (Kobayashi  and  Poff  1994)  is  modified  to  compute  wave  runup  on  the 
slope  which  is  assumed  to  be  impermeable.  The  details  of  the  computer  program 
VBREAK  have  been  published  separately  in  the  report  by  Kobayashi  and  Johnson 
(1995). 

In  Chapter  4,  the  model  VBREAK  is  initially  compared  with  a  previously 
developed  one- dimensional  model  IBREAK  (Kobayashi  and  Wurjanto  1989).  The 
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model  VBREAK  is  reduced  to  the  corresponding  one-dimensional  model  through 
the  appropriate  simplification  of  the  boundary  condition  for  the  momentum  flux 
correction.  The  two  models  are  then  based  on  the  same  governing  equations.  The 
MacCormack  method  is  assesed  using  one  of  the  tests  conducted  by  Ahrens  (1975) 
with  regular  waves  on  a  rough  1:2.5  slope.  Finally,  The  numerical  model  VBREAK 
is  compared  with  two  sets  of  regular  wave  data  in  Chapter  4.  One  data  set  is  the 
comprehensive  measurements  of  test  1  presented  by  Stive  (1980)  and  Stive  and  Wind 
(1982)  in  which  the  incident  regular  waves  broke  as  spilling  breakers  on  a  concrete 
1:40  beach.  The  other  data  set  is  the  detailed  velocity,  bottom  shear  stress  and  free 
surface  measurements  by  Cox  et  al.  (1995)  for  the  case  of  regular  waves  spilling  on 
a  rough,  impermeable  1:35  slope. 

The  summary  and  conclusions  of  this  report  are  given  in  Chapter  5 


Chapter  2 


MATHEMATICAL  FORMULATION 


2.1  Two-Dimensional  Equations  in  Shallow  Water 

The  approximate  governing  equations  adopted  in  the  numerical  model  VBREAK 
are  derived  from  the  two-dimensional  continuity  and  Reynolds  equations  (e.g.,  Rodi 
1980) 


dxj 

I  / 

df  ^  ^  dx\ 


=  0 


P  dx\  p  dx'j 


(2.1) 

(2.2) 


in  which  the  prime  indicates  the  physical  variables  and  the  summation  convention 
is  used  with  respect  to  repeated  indexes.  The  symbols  used  in  (2.1)  and  (2.2)  are 
depicted  in  Figure  2.1  where  t'  =  time;  =  horizontal  coordinate  taken  to  be  pos¬ 
itive  landward;  Xj  =  vertical  coordinate  taken  to  be  positive  upward  with  Xj  =  0  at 
the  still  water  level  (SWL);  u[  =  horizontal  velocity;  u'^  =  vertical  velocity;  p  =  fluid 
density  which  is  assumed  constant;  p'  =  pressure;  g  =  gravitational  acceleration;  6i2 
—  Kronecker  delta;  and  r-j  =  sum  of  turbulent  and  viscous  stresses.  Assuming  that 
the  viscous  stresses  are  negligible,  may  be  expressed  as  (e.g.,  Rodi  1980) 


R j  =  P 


Mx',  dx\]  3  “ 


(2.3) 


in  which  =  turbulent  eddy  viscosity;  and  k'  =  turbulent  kinetic  energy  per  unit 


Figure  2.1:  Definition  sketch. 


To  simplify  (2.1)  and  (2.2)  with  (2.3)  in  shallow  water,  the  dimensional  vari¬ 
ables  may  be  normalized  as 


T'  ’  ^  T'vW 

u[  U2 

^  vW  ’  ^  H'jT' 

v't  .  ,  _  k' 

H'yr  ’  y/^H'iT'  ’ 


a  =  r 


in  which  T'  and  H'  are  the  reference  wave  period  and  height  used  for  the  normal¬ 
ization.  The  parameter  a  defined  in  (2.6)  is  the  ratio  between  the  horizontal  and 
vertical  length  scales.  The  normalized  variables  in  (2.4)  and  (2.5)  are  assumed  to 


6 


be  on  the  order  of  unity  in  shallow  water.  The  normalization  of  v[  and  k'  in  (2.6)  is 
based  on  the  turbulence  measurements  in  a  wave  flume  by  Cox  et  al.  (1994)  which 
have  indicated  that  Ut  and  k  are  on  the  order  of  unity  or  less  inside  and  immediately 
outside  the  surf  zone,  respectively. 

Substituting  (2.4)-(2.6)  into  (2.1)-(2.3),  the  normalized  continuity  and  mo¬ 
mentum  equations  are  obtained.  The  conventional  notations  of  a;  =  a;i,  z  —  X2,  u  = 
u\  and  w  =  U2  are  used  in  the  following.  The  normalized  continuity  equation  is  given 

by 

du  dw 

The  momentum  equations  are  simplified  under  the  assumption  of  <7^  1  for  shallow 

water  waves.  The  approximate  horizontal  momentum  equation  is  expressed  as 


du  du  du  d  f  2k\  dr 


with 


T  =  Vt 


du 

dz 


The  approximate  vertical  momentum  equation  is  written  as 

d  (  2k\ 


(2.8) 


(2.9) 


(2.10) 


The  free  surface  and  bottom  are  located  at  z'  =  rj'  and  z'  =  as  shown  in 
Figure  2.1  where  the  bottom  is  assumed  to  be  fixed  and  impermeable.  The  water 
depth  h'  is  given  bj^  h'  =  {rj'  —  z'^).  The  dimensional  variables  rf',  z'^  and  h'  are 
normalized  by  the  vertical  length  scale  H' 


V  =  -TT, 


Zb  = 


H> 


(2.11) 


H'  '  ^  H' 

The  kinematic  boundary  conditions  at  the  free  surface  and  bottom  are  expressed  as 


dt]  dr) 

—  +  U-r - W  =  0  at  2  =  7/ 

dt  dx 

dzb 

u  — - u;  =  0  at  2  =  2f, 

dx 


(2.12) 

(2.13) 


The  normal  and  tangential  stresses  at  the  free  surface  are  assumed  to  be  zero.  These 
boundary  conditions  for  1  can  be  shown  to  yield 


0  at  z  =  Tj 


r  =  0  at  z  =  T] 


(2.14) 

(2.15) 


Integration  of  (2.10)  with  respect  to  z  using  (2.14)  gives 

2k 


p  =  7]  —  Z  — 
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(2.16) 


The  pressure  is  approximately  hydrostatic  in  shallow  water  where  k  is  on  the  order 
of  unity  or  less  and  a  is  relatively  large  to  satisfy  cr^  >  1.  Substituting  (2.16)  into 
(2.8),  the  horizontal  momentum  equation  is  rewritten  as 


du  du  du  dr]  dr 


(2.17) 


Eqs.  (2.7)  and  (2.17)  together  with  (2.9),  (2.12),  (2.13)  and  (2.15)  may  be  solved  nu¬ 
merically  but  considerable  numerical  difficulties  are  expected  because  the  unknown 
free  surface  elevation  r]  varies  rapidly  in  space  and  time.  In  addition,  such  a  numer¬ 
ical  model  will  be  too  time-consuming  to  compute  breaking  wave  motions  of  long 
duration. 


2.2  Depth-Integrated  Equations 

To  reduce  computational  efforts  significantly,  (2.7)  and  (2.17)  are  integrated 
from  z  =  Zb  to  z  =  T)  using  (2.12),  (2.13)  and  (2.15).  No  additional  approximation  is 
introduced  in  this  integration.  The  depth-integrated  continuity  equation  is  expressed 


as 


^+®i  =  o 

dt  dx 


(2.18) 


where  h  =  water  depth  given  by  h  =  {rj  —  Zb)’,  and  q  =  volume  flux  per  unit  width 


defined  as 


(2.19) 
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(2.20) 


The  depth-integrated  horizontal  momentum  equation  is  written  as 

with 

m=  f  [u  —  UY  dz  (2.21) 

Jzb 

in  which  U  =  depth- averaged  horizontal  velocity  defined  asU  =  q/h]  9  =  normalized 
bottom  slope  defined  as  0  =  dzbfdx',  u  =  bottom  shear  stress;  and  m  =  momentum 
flux  correction  due  to  the  vertical  variation  of  the  horizontal  velocity  u  where  m  =  0 
if  u  =  C/. 

The  previous  one-dimensional  models  IBREAK  (Kobayashi  and  Wurjanto 
1989),  RBREAK  (Wurjanto  and  Kobayashi  1991),  and  RBREAK2  (Kobayashi 
and  Poff  1994)  assumed  m  =  0  and  expressed  tj  in  terms  of  U .  Eqs.  (2.18)  and 
(2.20)  with  m  =  0  were  solved  using  the  dissipative  Lax-Wendroff  finite  difference 
method  to  compute  h  and  9  as  a  function  of  t  and  x.  These  one-dimensional  models 
do  not  predict  the  vertical  variations  of  the  fluid  velocities  u  and  w.  Furthermore, 
these  models  do  not  account  for  energy  dissipation  due  to  wave  breaking  explicitly. 

Boussinesq  equations  have  been  extended  to  predict  breaking  waves  on  gentle 
slopes  (Zelt  1991;  Schaffer  et  al.  1992).  Boussinesq  equations  without  the  dispersive 
terms  correspond  to  (2.18)  and  (2.20)  if  the  bottom  friction  is  included  in  Boussi¬ 
nesq  equations  (Zelt  1991).  Gharangik  and  Chaudhry  (1991)  computed  hydraulic 
jumps  using  Boussinesq  equations  with  and  without  the  dispersive  terms  and  found 
that  the  dispersive  terms  had  little  effect  on  the  computed  hydraulic  jumps.  This 
indicates  that  the  dispersive  terms  may  be  negligible  for  breaking  waves  inside  the 
surf  zone.  Moreover,  the  dispersive  terms  derived  under  the  assumption  of  potential 
flow  may  not  be  valid  for  breaking  waves.  To  include  energy  dissipation  due  to  wave 
breaking  in  Boussinesq  equations,  Zelt  (1991)  and  Schaffer  et  al.  (1992)  added  a 
term  corresponding  to  the  term  for  the  momentum  flux  correction  m  in  (2.20).  Zelt 
(1991)  expressed  this  additional  term  in  the  form  of  horizontal  momentum  diffusion 


9 


with  an  artificial  viscosity  proposed  by  Heitner  and  Housner  (1970).  The  artificial 
viscosity  was  calibrated  for  breaking  solitary  waves  where  the  diffusion  term  was 
activated  using  a  semi-empirical  criterion  for  solitary  wave  breaking.  On  the  other 
hand,  Schaffer  et  al.  (1992)  expressed  the  additional  momentum  flux  using  a  simple 
approach  based  on  a  surface  roller  that  represented  a  passive  bulk  of  water  riding  on 
the  front  of  a  breaking  wave.  An  empirical  geometric  method  was  used  to  determine 
the  shape  and  location  of  the  surface  rollers  during  the  computation.  These  models 
do  not  predict  the  A^ertical  variations  of  the  fluid  velocities.  It  is  also  not  certain 
whether  the  computed  energy  dissipation  was  truly  caused  by  the  term  added  to 
the  momentum  equation  because  they  did  not  check  whether  the  computed  results 
satisfied  the  energy  equation  as  will  be  elaborated  in  the  following. 

In  this  report,  the  equation  for  the  momentum  flux  correction  m  is  derived 
from  the  depth-integrated  instantaneous  wave  energy  equation  (Kobayashi  and  Wur- 
janto  1992) 

f  +  =  (2,22) 

which  is  obtained  by  integrating  (2.17)  multiplied  by  u  from  z  =  Zb  to  z  =  rj  hy 
use  of  (2.12),  (2.13),  (2.15)  and  (2.18).  The  specific  energy  E  defined  as  the  sum  of 
kinetic  and  potential  energy  per  unit  horizontal  area  is  given  by 

E  =  ^{qU  +  m  +  rj^)  for  zj,  <  0  (2.23) 

^  ~  \  ^  -  -s^t)  for  Zb>  0  (2.24) 

in  which  the  potential  energy  is  taken  to  be  relative  to  the  potential  energy  in  the 

absence  of  wave  action  with  SWL  at  =  0.  The  energy  flux  Ep  per  unit  width  is 

expressed  as 

Ep  =  r}q  +  ^  (^qU^  -(-  3mU  -|-  m3)  (2.25) 

with 

^3=  f  {u  —  U)^  dz  (2.26) 

Et 
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in  which  m3  =  kinetic  energy  flux  correction  due  to  the  third  moment  of  the  velocity 
deviation  (u  —  U)  over  the  depth  where  m3  =  0  if  u  =  t/.  The  energy  dissipation 
rate  D  per  unit  horizontal  area  in  (2.22)  is  given  by 

D  =  [\^dz  (2.27) 

Jzb  dz 

where  use  is  made  of  the  no  slip  condition  u  =  0  a.t  z  =  Zb- 

The  wave  boundary  layer  is  not  analyzed  explicitly  in  this  numerical  model. 
The  energy  dissipation  rate  Df  inside  the  wave  boundary  layer  may  be  estimated 
by  (Jonsson  and  Carlsen  1976) 

Df  =  TbUb  (2.28) 

where  Ub  =  near-bottom  horizontal  velocity  immediately  outside  the  wave  boundary 
layer.  The  normalized  bottom  shear  stress  T(,  may  be  expressed  as 

n  =  fw\ub\ub  ;  ^  ^f'w  (2-29) 

in  which  fl,  =  wave  friction  factor  (Jonsson  1966).  The  value  of  f!^  specified  as 
input  is  allowed  to  vary  spatially  to  accommodate  the  spatial  variation  of  bottom 
roughness  (Kobayashi  and  Raichle  1994).  The  previous  one- dimensional  models 
(e.g.,  Kobayashi  and  Wurjanto  1992)  employed  (2.27)  and  (2.28)  in  which  the  depth- 
averaged  velocity  U  and  the  corresponding  friction  factor  /'  were  used  instead  of  the 
near-bottom  velocity  Ub  and  the  wave  friction  factor  /(^.  Cox  et  al.  (1995)  showed 
that  the  bottom  shear  stress  and  near-bottom  velocity  measured  inside  the  surf  zone 
could  be  related  fairly  well  by  the  quadratic  friction  equation  (2.28)  with  the  wave 
friction  factor  estimated  using  the  formula  of  Jonsson  (1966). 

The  energy  dissipation  rate  D  given  by  (2.26)  may  be  expressed  as 


D  =  Df  J-  Db 


(2.30) 


in  which  Db  =  energy  dissipation  rate  outside  the  wave  boundary  layer  due  to  wave 
breaking.  Assuming  that  the  thickness  of  the  wave  boundary  layer  is  much  smaller 
than  the  water  depth,  Db  may  be  estimated  using  (2.26)  together  with  (2.9) 


Db  = 


outside  boundary  layer 


(2.31) 


where  the  vertical  variations  of  u  and  i/t  outside  the  wave  boundary  layer  will  be 
assumed  in  the  following. 

Rearranging  the  instantaneous  wave  energy  equation  (2.22)  with  (2.27)  and 
(2.29)  by  use  of  (2.18)  and  (2.20),  the  equation  for  the  momentum  flux  correction 
m  is  derived 


dm  d 


din 

2U  — - 2{TbUb  +  Db) 

ox 


(2.32) 


with 


Ub  =  Ub  —  U 


(2.33) 


in  which  Ub  =  near-bottom  horizontal  velocity  correction  due  to  the  vertical  variation 
of  the  horizontal  velocity  u  outside  the  wave  boundary  layer.  If  u  =  f/,  {tj  =  0,  m  =  0 
and  m3  =  0.  Asa  result,  (2.31)  yields  Db  =  0  ifu  =  U,  whereas  Db  given  by  (2.30) 
outside  the  wave  boundary  layer  is  also  zero  if  u  is  independent  of  z.  This  proves 
that  the  energy  dissipation  due  to  wave  breaking  in  the  previous  one-dimensional 
models  based  on  the  assumptions  of  W6  =  0,  m  =  0  and  m3  =  0  is  solely  numerical 
(Kobayashi  and  Wurjanto  1992). 

In  order  to  express  m,  m3  and  Db  in  terms  of  Ub,  the  horizontal  velocity  u 
outside  the  wave  boundary  layer  is  assumed  to  be  expressible  in  the  form 


u{t,  x,z)  =  U (t,  x)  -t-  Ub(t,  x)F(C)  (2.34) 

with 


C  =  [z  —  2)6(x)]  /h{t,  x)  for  0  <  C  <  1 


(2.35) 


in  which  F  =  normalized  function  expressing  the  vertical  variation  of  the  velocity- 
deviation  (u  —  U)  from  (  =  0  immediately  outside  the  wave  boundary  layer  to  C  =  1 
at  the  free  surface.  Furthermore,  the  dimensional  turbulent  eddy  viscosity  v[  outside 
the  wave  boundary  layer  is  assumed  to  be  given  by 

<  =  (C,hf  ^  (2.36) 

in  which  Ce  =  mixing  length  parameter.  The  turbulence  measurements  inside  the 
surf  zone  by  Cox  et  al.  (1994)  have  indicated  that  (2.35)  is  a  reasonable  first  approx¬ 
imation  outside  the  wave  boundary  layer  and  that  Ci  is  on  the  order  of  0.1.  Using 
(2.4)-(2.6)  and  (2.11),  the  normalized  turbulent  eddy  viscosity  vt  corresponding  to 
(2.35)  is  expressed  as 

vt  =  Cja  (2.37) 

Substitution  of  (2.33)  with  (2.34)  and  (2.36)  into  (2.21),  (2.25)  and  (2.30)  yields 


rn  =  C2hul 

;  C2  = 

f 

(2.38) 

1713  =  Czhul 

;  C3  = 

F^dC 

Jo 

(2.39) 

F  dF  ^ 

Db  =  CBCla\uif  : 

;  Cb  — 

Jo  dC 

(2.40) 

in  which  m  and  Db  are  positive  or  zero. 

Madsen  and  Svendsen  (1983)  and  Svendsen  and  Madsen  (1984)  assumed  a 
cubic  velocity  profile  for  their  analyses  of  a  hydraulic  jump  and  a  turbulent  bore  on 
a  beach.  Accordingly,  the  function  F  in  (2.33)  outside  the  wave  boundary  layer  is 
assumed  to  be  cubic  and  expressed  as 

T’  =  1  —  (3  +  0.75a)(^^  -f  for  0  <  C  <  1  (2.41) 

in  which  a  —  cubic  velocity  profile  parameter.  The  function  F  given  by  (2.40) 
satisfies  (2.19)  with  q  =  Uh  and  (2.32).  The  shear  stress  t  given  by  (2.9)  with 
(2.36)  must  satisfy  (2.15).  However,  (2.40)  yields  t  =  0  at  C  =  1  only  if  a  =  4. 
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Moreover,  (2.40)  results  in  r  =  0  at  (^  =  0  immediately  outside  the  wave  boundary 
layer  in  contradiction  with  the  turbulence  measurements  inside  the  surf  zone  by 
Cox  et  al.  (1994).  Consequently,  (2.40)  with  the  single  empirical  parameter  a 
may  not  predict  the  shear  stress  accurately  in  the  vicinity  of  the  free  surface  and 
bottom.  Comparison  of  (2.40)  with  the  cubic  profile  assumed  by  Svendsen  and 
Madsen  (1984)  suggests  that  the  parameter  a  is  approximately  3.  The  range  of  a 
=  3-4  is  considered  in  the  following.  Substitution  of  (2.40)  into  the  equations  for 
C2,  Cz  and  Cb  in  (2.37)-(2.39)  yields 


^2 

Cz 

Cb 


2b  a  IP'  ah  aP 

3a  36^  3a^  +  daP  aPb  aP 


(2.42) 

(2.43) 

(2.44) 


in  which  6  =  —  (3  +  0.75a). 

Figure  2.2  shows  the  cubic  velocity  profile  function  F  given  by  (2.40)  as  a 
function  of  ^  for  a  =  3.0,  3.5  and  4.0.  The  abscissa  in  Figure  2.2  is  the  value  of 
—F  because  wj  in  (2.33)  is  expected  to  be  negative  under  the  wave  crest.  Figure  2.2 
hence  depicts  the  normalized  vertical  variation  of  the  horizontal  velocity  deviation 
{u  —  U)  under  the  wave  crest.  The  assumed  cubic  profile  is  not  sensitive  to  the 
parameter  a  in  the  range  of  a  =  3-4  except  in  the  vicinity  of  the  free  surface  where 
no  velocity  data  is  available  inside  the  surf  zone.  Figure  2.3  shows  the  parameters 
C2,  Cz  and  Cb  as  a  function  of  the  cubic  profile  parameter  a.  These  parameters  vary 
little  for  a  =  3-4.  Figure  2.3  indicates  that  C2  —  0.5,  Cz  —  —0.03  and  Cb  —  13. 
In  short.  Figures  2.2  and  2.3  imply  that  the  computed  results  will  not  be  sensitive 
to  the  empirical  parameter  a.  The  mixing  length  parameter  C(  affects  only  Db 
given  by  (2.39)  but  will  modify  the  computed  magnitude  of  Db  more  than  the  cubic 
profile  parameter  a  because  Db  is  proportional  to  C^. 
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Figure  2.3:  C21  Cz,  and  Cb  as  a  function  of  cubic  profile  parameter 


Eqs.  (2.18),  (2.20),  and  (2.31)  together  with  (2.28),  (2.32),  and  (2.37)-(2.39) 
will  be  solved  numerically  in  the  next  section  to  compute  h,  q  and  m  as  a  function 
of  t  and  X.  To  obtain  Ub  using  (2.37)  for  the  computed  h  and  m,  it  is  assumed  that 


/mV 

^6  = 

\C2h) 

/  m  y/2 

Ub  = 

Kc^h) 

for  17  >  0 
for  17  <  0 


(2.45) 

(2.46) 


which  ensures  that  |u6|  <  |17|  with  Ub  =  (17  +  ut).  It  is  required  in  (2.44)  that  m  >  0. 
For  the  computed  h,  U  =  qjh  and  ut,  the  horizontal  velocity  u  can  be  obtained 
using  (2.33).  The  vertical  velocity  w  can  be  found  using  the  continuity  equation 

in  which  (  is  given  by  (2.34)  and  use  is  made  of  w  =  0  z  =  Zb- 

To  examine  the  degree  of  numerical  dissipation  hidden  in  the  computed  re¬ 
sults,  the  instantaneous  energy  equation  (2.22)  with  (2.29)  is  averaged  from  t  =  <stat 
fot  =  linax 


AE  (^Epj  =  —Df  —  Db 


with 


AE 


dx 


E{t  =  tiaax)  -  E{t  =  tsut) 


(2.48) 

(2.49) 


Imax  Istat 

in  which  the  overbar  denotes  the  time  averaging  from  the  starting  time  tsta.t  of 
the  statistical  calculations  to  the  ending  time  tmax  of  the  computation  as  will  be 
explained  in  Section  3.4.  For  the  computed  h,  q  and  m,  E,  Ef,  Df  and  Db  are 
computed  using  (2.23),  (2.24),  (2.27)  and  (2.39),  respectively,  during  igtat  <  tmax- 
The  computed  E,  Ep,  Df  and  Db  will  satisfy  the  time-  averaged  energy  equation 
(2.46)  in  the  absence  of  numerical  dissipation  in  the  adopted  numerical  procedures. 
In  the  previous  one- dimensional  models,  Db  was  calculated  using  (2.46)  because 
these  models  did  not  include  any  physical  dissipation  mechanism  associated  with 
wave  breaking. 


Chapter  3 


NUMERICAL  MODEL 


3.1  MacCormack  Method 

To  solve  (2.18),  (2.20)  and  (2.31)  for  /«,  q  and  m,  these  equations  are  com¬ 
bined  and  expressed  in  the  following  vector  form: 


with 


and 


dV  dF  ^  ^ 


h 

0 

9 

;  F  = 

F2 

Q 

II 

G2 

m 

1 

1 _ 

1 

CO 

1 _ 

F2 

Fs 


qU  +  m  + 
3mU  -|-  m3 


G2  =  Oh  Tf) 


;  G3  =  2  (  TbUb  +  Db  —  U 


dm' 
dx  , 


(3.1) 


(3.2) 


(3.3) 

(3.4) 


in  which  U  —  qjh.  rj  is  given  by  (2.28)  with  =  (17  +  Ub)  and  Ub  is  calculated 
using  (2.44).  m3  and  Db  are  obtained  from  (2.38)  and  (2.39),  respectively.  Eq. 
(  3.1)  is  solved  numerically  using  the  MacCormack  method  (MacCormack  1969) 
which  is  a  simplified  variation  of  the  two-step  Lax-Wendroff  method  (e.g.,  Anderson 
et  al.  1984)  and  has  been  applied  successfully  for  the  computation  of  unsteady  open 
channel  flows  with  hydraulic  jumps  (e.g.,  Fennema  and  Chaudhry  1986;  Gharangik 


and  Chaudhry  1991).  The  use  of  the  Lax-Wendroff  method  for  (  3.1)  would  be  very 
difficult  because  this  method  requires  the  Jacobian  of  F  with  respect  to  U. 

The  initial  time  t  =  0  for  the  computation  marching  forward  in  time  is  taken 
to  be  the  time  when  the  incident  wave  arrives  at  the  seaward  boundary  located  at 
X  =  0  and  there  is  no  wave  action  in  the  computation  domain  a:  >  0.  The  initial 
conditions  for  the  computation  are  thus  given  by 


at  t  =  0  for  Z(,  <  0 

(below  SWL) 

(3.5) 

at  t  =  0  for  Zfr  >  0 

(above  SWL) 

(3.6) 

=  0  ;  m  =  0  at  t 

=  0 

(3.7) 

in  which  Zb  =  normalized  bottom  elevation  taken  to  be  negative  below  SWL. 

A  finite  difference  grid  of  constant  nodal  interval  Ax  and  variable  time  step 
At  is  used  in  the  numerical  model  VBREAK.  The  spatial  nodes  are  located  at 
X  =  {j  —  l)Ax  with  j  =  I,  2,  . . . ,  jmax  where  jmax  =  number  of  the  spatial  nodes 
in  the  computation  domain.  The  computational  shoreline  is  defined  as  the  location 
where  the  normalized  instantaneous  water  depth  h  equals  a  small  value  6  such  as 
S  =  10“^  as  in  the  previous  one-dimensional  models  (e.g.,  Kobayashi  and  Poff  1994). 
The  integer  s  is  used  to  indicate  the  wet  node  next  to  the  moving  shoreline  such  that 
hs+i  ^  ^  <  hs  where  hs  and  hg+i  are  the  values  of  hj  at  the  node  j  =  s  and  (s  -f-l), 
respectively.  It  is  noted  that  no  wave  overtopping  and  transmission  are  allowed  in 
the  present  form  of  VBREAK  unlike  the  previous  one-dimensional  models.  It  is 
hence  required  that  s  <  jmax- 

The  values  of  Uj  at  the  node  j  with  j  =  1,  2,  . . . ,  s  and  at  the  present  time 
t  are  known  in  the  following  where  Uj  =  0  with  j  =  (s  +  1),  (s  -|-  2),  . . . ,  jmax 
landward  of  the  shoreline  node  s.  The  unknown  values  of  U!  at  the  node  j  and 

J  *' 


at  the  next  time  level  t*  =  {t  +  At)  are  denoted  by  the  superscript  asterisk.  The 
predictor,  corrector  and  final  steps  of  the  MacCormack  method  are  expressed  as 


Ui 

=  u,  -  ^  (Frt.  -  F,-)  -  A(G, 

for  j  =  1,  2,  . . . ,  s 

(3.8) 

Ui 

=  U,-^(f,-F,--i)-A(G, 

for  =  2,  3,  . . . ,  s 

(3.9) 

u* 

=  5(u,  +  u,) 

for  j  =  2,  3,  . . . ,  s 

(3.10) 

in  which  a  forward  spatial  difference  is  used  for  the  second  term  on  the  right  hand 
side  of  the  predictor  equation  (  3.8),  and  a  backward  spatial  difference  is  used  for 
this  term  in  the  corrector  equation  (  3.9).  This  is  because  the  spatial  difference 
in  the  predictor  equation  is  recommended  to  be  in  the  direction  of  propagation  of 
wave  fronts  (Anderson  et  al.  1984).  Accordingly,  the  term  dmidx  in  the  function 
Gs  defined  in  (  3.4)  is  expressed  by  the  forward  and  backward  spatial  differences 
in  (  3.8)  and  (  3.9),  respectively.  The  values  of  Uj  computed  by  (  3.8)  and  the 
corresponding  values  of  Fj  and  Gj  in  (  3.9)  are  temporary  values  at  the  next  time 
level  t*.  In  the  computer  program  VBREAK,  the  three  equations  corresponding  to 
each  of  (  3.8),  (  3.9)  and  (  3.10)  are  used  for  clarity.  The  values  of  are  computed 
using  the  seaward  boundary  conditions  in  Section  3.5.  The  numerical  procedures 
for  the  moving  shoreline  in  Section  3.6  are  used  to  improve  the  computed  values  of 
U*  and  find  the  shoreline  node  s*  at  the  next  time  level  t*. 

3.2  Numerical  Stability  and  Smoothing 

The  constant  nodal  interval  Ax  needs  to  be  small  enough  to  resolve  steep 
wave  fronts  in  the  siirf  zone.  The  variable  time  step  size  At  for  numerically  stable 
computation  is  estimated  at  the  beginning  of  each  time  step  using  the  following 
approximate  equation: 

C  At 

" - /irri  ,  fr\  for  i  =  1,  2,  . . . ,  s  (3.11) 

max(|C/j|  +  Jhj) 
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in  which  Cn  is  the  Courant  number  and  the  denominator  in  (  3.11)  is  the  maximum 
value  of  (|17,  |  +  at  all  the  wet  nodes  at  the  present  time  t.  The  value  of  At 
for  each  time  step  is  selected  using  (  3.11)  except  for  the  last  time  step  that  is 
chosen  such  that  t*  —  {t  +  At)  =  t^iax  for  given  t  and  tmax-  The  numerical  stability 
of  the  MacCormack  method  applied  to  (  3.1)  with  m  =  0  and  G2  =  0  requires 
that  C'n  <  1  (e.g.,  Anderson  et  al.  1984).  Eq.  (  3.11)  is  approximate  because  the 
characteristic  equations  corresponding  to  (  3.1)  with  m  7^  0  can  not  be  expressed 
in  simple  analytical  forms.  Moreover,  (  3.11)  does  not  account  for  the  shoreline 
algorithm  which  tends  to  suffer  numerical  difficulties.  Consequently,  the  value  of 
Cn  less  than  unity  is  specified  as  input  to  adjust  At  for  the  successful  computation 
of  each  case.  This  manual  adjustment  of  At  through  the  specification  of  Cn  on  the 
order  of  0.4  appears  to  be  sufficient  because  the  computation  time  of  VBREAK  is 
relatively  short. 

Use  of  the  MacCormack  method  results  in  numerical  high-frequency  oscilla¬ 
tions  which  tend  to  appear  at  the  rear  of  a  breaking  wave,  especially  on  a  gentle 
slope.  For  open-channel  flows,  Chaudhry  (1993)  summarized  a  procedure  presented 
by  Jameson  et  al.  (1981)  to  smooth  these  high-frequency  oscillations  without  dis¬ 
turbing  the  rest  of  the  computed  variations.  To  apply  this  procedure  for  breaking 
waves  on  slopes,  the  computed  water  depth  at  the  node  j  and  at  the  next  time 
level  t*  is  used  to  calculate  the  parameter  Vj  at  the  node  j  defined  as 

\h’^,-2h-  +  hU\ 

’  l'‘•+ll+2|/^'|  +  |/!•_,|  ^  ’  . *  *  * 

The  parameter  Cj+o.s  at  the  midpoint  of  the  nodes  j  and  (i  +  1)  is  given  by 


(3.12) 


K + I'h 


max(z/j,  I'j+i)  for  j  =  2,  3,  ...,  (s*  —  2)  (3.13) 


in  which  k  =  numerical  damping  coefficient  for  regulating  the  amount  of  damping 
the  high-frequency  oscillations.  The  computed  water  depth  h*-  is  modified  as 

=  /j^  +  ej+0.5  -  /ij)  -ej-0.5  (h*j  -  ^*_i)  for  j  =  3,  4,  . . . ,  (s*-2)  (3.14) 
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which  should  be  considered  as  a  FORTRAN  replacement  statement.  Likewise,  UJ 
and  are  smoothed  using  (  3.14)  with  h’j  being  replaced  by  C/J  and  m*,  respectively, 
where  Cj+o.s  is  the  same.  The  smoothed  h’j  and  Uj  are  used  to  calculate  q'j  =  hjUj. 

Chaudhry  (1993)  suggested  the  expressions  of  vj  at  the  end  points  j  =  I 
and  s*  in  (  3.12)  for  open-channel  flows.  Addition  of  these  expressions  in  (  3.13) 
and  (  3.14)  is  found  to  produce  spurious  fluid  motions  even  in  the  absence  of  waves 
on  slopes.  As  a  result,  the  smoothing  at  the  end  points  is  not  recommended  for 
breaking  waves  on  slopes. 

The  numerical  damping  coefficient  k  is  specified  as  input  to  VBREAK.  For 
breaking  waves  on  gentle  beach  slopes,  the  value  of  k  on  the  order  of  unity  is  found  to 
be  necessary  to  damp  the  high-frequency  oscillations  adequately.  For  waves  surging 
on  steep  slopes  of  coastal  structures,  the  value  of  k  on  the  order  of  0.1  appears  to 
be  sufficient.  However,  the  smoothing  procedure  based  on  (  3.12)  tends  to  cause 
more  damping  near  the  shoreline  where  the  water  depth  h  is  very  small.  To  remedy 
this  uneven  damping,  the  term  [{h’j  -f  /ij+i)/2]°-®  is  added  in  (  3.13)  to  reduce  the 
damping  near  the  shoreline. 

3.3  Incident  or  Measured  Wave  Profile 

The  required  wave  input  for  the  numerical  model  VBREAK  is  the  normalized 

incident  or  measured  wave  profile  at  the  seaward  boundary  of  the  computation 

domain,  that  is,  Tji{t)  =  r]l{t')fH'  or  T]{t)  =  T]'{t')fH'  at  a;  =  0  with  t  =  t'/T'  where 

H'  and  T'  are  the  reference  wave  height  and  period  used  for  the  normalization  in 

(2.4)-(2.6).  The  specification  of  T]i{t)  or  j/(f)  at  a;  =  0  for  t  >  0  needs  to  satisfy  the 

condition  ?/,•  =  0  or  ?/  =  0  at  a:  =  0  at  the  initial  time  f  =  0  to  be  consistent  with 
« 

the  assumed  initial  conditions  of  no  wave  action  in  the  region  of  a;  >  0  at  t  =  0. 

The  temporal  variation  of  the  incident  wave  profile  T}i{t)  at  a;  =  0  can  be 
1.  the  incident  monochromatic  wave  profile  computed  (by  the  computer  program 
VBREAK)  using  an  appropriate  wave  theory,  or 
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2.  a  user-specified  incident  irregular  or  transient  wave  train  including  the  inci¬ 
dent  wave  profile  measured  in  the  absence  of  wave  reflection,  measured  in  the 
presence  of  wave  reflection  but  separated  from  the  reflected  wave  using  linear 
wave  theory,  or  generated  numerically  for  given  frequency  spectrum. 

For  convenience,  the  former  will  be  referred  to  as  the  case  of  regular  waves,  while 
the  latter  is  simply  called  the  case  of  irregular  waves,  although  this  case  can  also 
include  monochromatic  transient  waves. 

The  temporal  variation  of  the  normalized  free  surface  elevation  ?/(t)  at  x  =  0 
can  be  specified  as  input  if  r}{t)  at  a;  =  0  is  measured  in  the  presence  of  wave  reflec¬ 
tion.  This  option  eliminates  uncertainties  associated  with  the  separation  of  incident 
and  reflected  waves  using  linear  wave  theory  in  laboratory  and  field  measurements. 

3.3.1  Incident  Regular  Wave 

For  the  case  of  incident  regular  waves,  the  periodic  variation  of  is  computed  by 
the  computer  program  using  either  cnoidal  or  Stokes  second-order  wave  theory.  The 
height  and  period  of  the  incident  regular  waves  at  the  seaward  boundary  located  at 
a;  =  0  are  denoted  by  HI  and  T^.  The  reference  wave  period  T'  is  taken  as  T'  =  T/ 
for  the  incident  regular  waves.  The  reference  wave  height  H'  specified  as  input  may 
be  in  deeper  water  depth.  Since  the  numerical  model  is  based  on  the  assumption  of 
shallow  water  waves,  the  seaward  boundary  should  be  located  in  relatively  shallow 
water.  As  a  result,  it  is  not  always  possible  to  take  H'  =  H'-.  Defining  Ks  =  H'JH', 
the  height  and  period  of  the  regular  wave  profile  at  a:  =  0  is  ilTs  and  unity, 
respectively. 

For  Stokes  second-order  wave  theory,  the  incident  wave  profile  T]i{t)  at  a;  =  0 
is  given  by  (e.g.,  Kobayashi  and  Poff  1994) 

Tji{t)  =  Ks  I  +  to)]  +  a2Cos  -1-  to)]|  for  t>0  (3.15) 
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L  —  Lq  tank—  (3-17) 

Jj 

where  to  =  time  shift  computed  to  satisfy  the  conditions  that  ?/,■  =  0  at  t  =  0 
and  T]i  decreases  initially;  02  =  normalized  amplitude  of  the  second-order  harmonic; 
L  =  L'ld[]  dt  =  d^lH'\  Lq  =  L'oM't'i  —  water  depth  below  SWL  at  a:  =  0; 
L'  =  dimensional  wavelength  at  a;  =  0;  and  Lq  =  dimensional  linear  wavelength 
in  deep  water.  The  normalized  wavelength  L  satisfying  (  3.17)  for  given  Lq  is 
computed  using  a  Newton-Raphson  iteration  method.  Eq.  (  3.16)  yields  the  value 
of  <i2  for  given  dt  =  d't/H',  Kg.,  and  L.  Since  (  3.15)  satisfies  T)i{t  -|-  1)  =  T]i{t) 
and  —  to)  =  r)i{t  to),  it  is  sufficient  to  compute  the  profile  t]i{t)  for  0  < 

{t  +  to)  <  0.5  to  obtain  rji{t)  for  0  <  t  <  tmax  where  tmax  =  specified  computation 
duration.  Eq.  (  3.15)  may  be  appropriate  if  the  Ursell  parameter  Ur  <  26  where 
Ur  =  [Hl(L')‘^ / [d't)^]  =  jdt)  at  a;  =  0.  It  is  noted  that  the  value  of  Ur  based 

on  the  normalized  wavelength  L  computed  from  (  3.17)  is  simply  used  to  decide 
whether  cnoidal  or  Stokes  second-order  wave  theory  is  applied. 

For  the  case  of  Ur  >26,  cnoidal  wave  theory  is  used  to  compute  the  incident 
wave  profile  'qi{t)  at  x  =  0  (e.g.,  Kobayashi  and  Poff  1994) 


»?i(0  =  'Hmin  +  Ks  cn^  [2if(t  -|-  to)]  for  f  >  0  (3.18) 

with 

Vmin  =  “  (1  -  f-)  -  (3-19) 

where  rjiain  —  normalized  trough  elevation  below  SWL;  cn  =  Jacobian  elliptic  func¬ 
tion;  K  =  complete  elliptic  integral  of  the  first  kind;  E  =  complete  elliptic  integral 
of  the  second  kind,  which  should  be  differentiated  from  the  specific  wave  energy 
E  given  by  (2.23);  and  m  =  parameter  determining  the  complete  elliptic  integrals 
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K{m)  and  E{m).  It  is  noted  that  this  parameter  m  is  different  from  the  momen¬ 
tum  flux  correction  m  used  in  the  previous  sections.  The  notation  of  m  for  cnoidal 
wave  theory  is  standard  and  therefore  used  here.  The  parameter  m  for  cnoidal  wave 
theory  is  related  to  the  Ursell  parameter  Ur 

Ur  =  (3.20) 

dt  3 


For  Ur  >  26,  the  parameter  m  is  in  the  range  0.8  <  m  <  1.  The  parameter  m  for 
given  (T,  dt,  and  Kg,  where  a  is  defined  in  (2.6),  is  computed  from 


1  + 


Ks 
m  dt 


m  +  2  — 


-1  =  0 


(3.21) 


where  the  normalized  wavelength  L  is  given  by  (  3.20)  as  a  function  of  m  for  given 
dt  and  Kg.  The  left  hand  side  of  (  3.21)  is  a  reasonably  simple  function  of  m  in  the 
range  0.8  <  m  <  1.  As  a  result,  (  3.21)  can  be  solved  using  an  iteration  method 
which  successively  narrows  down  the  range  of  m  bracketing  the  root  of  (  3.21). 
After  the  value  of  m  is  computed  for  given  a,  dt,  and  Kg,  the  values  of  Ur  and  L 
are  computed  using  (  3.20),  while  (  3.19)  yields  the  value  of  rjmin-  The  incident 
wave  profile  T)i{t)  is  computed  using  (  3.18)  for  0  <  (t  +  io)  <  0.5  where  the  time 
shift  to  and  the  periodicity  and  symmetry  of  the  cnoidal  wave  profile  are  used  in  the 
same  manner  as  the  Stokes  second-order  wave  profile  given  by  (  3.15).  It  should  be 
mentioned  that  the  Jacobian  elliptic  function  and  the  complete  elliptic  integrals  of 
the  first  and  second  kinds  are  computed  using  the  subroutines  given  by  Press  et  al. 
(1986). 


3.3.2  Incident  Irregular  Wave 

The  specification  of  incident  irregular  waves  as  input  to  VBREAK  is  the  same 
as  in  the  previous  one-dimensional  model  RBREAK2  (Kobayashi  and  Poff  1994). 
Various  examples  of  user-specified  irregular  wave  trains  were  given  by  Wurjanto  and 
Kobayashi  (1991). 
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The  incident  wave  train  rii{t)  normalized  by  the  reference  wave  height  H'  is 
read  at  the  following  sampling  rate 

c,  _  ^max 

'  “  NPINP-1 

in  which  NPINP  =  specified  number  of  points  in  the  input  wave  train  for  0  <  f  <  fmax 
with  t  and  tmax  being  normalized  by  the  reference  wave  period  T'.  The  reference 
wave  height  and  period  can  be  chosen  as  any  height  and  period  that  are  convenient 
for  the  analysis  of  computed  normalized  results.  The  normalized  wave  height  Kg 
specified  as  input  for  case  of  incident  regular  waves  is  not  required  when  an  irregular 
wave  time  series  is  used.  It  is  noted  that  the  normalized  incident  regular  wave  train 
given  by  (  3.15)  and  (  3.18)  is  also  computed  at  the  rate  Sti  given  by  (  3.22). 

The  sampling  rate  6ti  must  be  small  enough  to  resolve  the  temporal  variation 
of  ?/t(t)  but  is  normally  much  larger  than  the  finite  difference  time  step  At  calculated 
by  (  3.11)  for  the  numerically  stable  computation.  A  simple  linear  interpolation  of 
T]i{t)  sampled  at  the  rate  6ti  is  performed  to  find  the  value  of  at  the  time  level 

t*  =  (t  +  At)  during  the  time-marching  computation. 

3.3.3  Measured  Wave  Profile 

If  the  free  surface  oscillation  is  measured  at  the  seaward  boundary  of  the 
computation,  it  is  more  direct  and  straightforward  to  specify  the  measured  free 
surface  oscillation  as  input  to  VBREAK  .  This  option  eliminates  the  uncertainty 
associated  with  the  separation  of  incident  and  reflected  waves  using  linear  wave 
theory. 

The  measured  free  surface  elevation  above  SWL  at  the  seaward  boundary  is 
normalized  by  the  values  of  H'  and  T'  specified  as  input.  The  normalized  input 
time  series  of  7/(i)  at  a;  =  0  is  read  at  the  sampling  rate  8ti  given  by  (  3.22)  and 
interpolated  linearly  in  the  same  way  as  7/,(t)  at  x  =  0  for  case  of  specified  irregular 


waves. 


3.4  Statistical  C^alculations 


The  statistical  calculations  in  this  report  imply  the  calculations  of  the  mean, 
root-mean-square  (rms),  maximum  and  minimum  values  of  computed  time-series  for 
the  duration  of  4tat  ^t  <  tmax-  For  example,  the  rms  value  of  77  is  defined  as 


^rms 


iv-v) 


1/2 


7/2  -  (nf 


1/2 


(3.23) 


in  which  the  overbar  denotes  the  time  averaging  for  4tat  <  ^  <  ^max- 

For  the  case  of  regular  waves,  the  statistical  calculations  should  be  executed 
over  the  last  wave  period,  assuming  that  the  computation  duration  tj^ax  is  large 
enough  to  reach  the  periodicity  of  time-varying  quantities  during  4tat  <  t  <  tmax 
with  tstat  =  (tmax  ~  !)•  This  duration  ranges  from  approximately  five  wave  periods 
for  coeistal  structures  (Kobayashi  and  Wurjanto  1989)  to  about  30  wave  periods  for 
beaches  (Kobayashi  et  al.  1989). 

For  irregular  wave  computations,  the  statistical  calculations  are  conducted 
over  most  or  all  of  the  computation  duration.  The  initial  transient  waves  may  be 
excluded  by  specifying  an  appropriate  value  of  istat  estimated  from  the  corresponding 
regular  wave  computation. 


3.5  Wave  Reflection 

The  normalized  reflected  wave  train  gr{t)  at  the  seaward  boundary  needs  to 
be  computed  to  estimate  the  degree  of  wave  reflection  from  the  computation  domain. 
It  is  also  required  to  find  the  unknown  value  of  the  vector  Uj  at  a;  =  0  and  at  the 
next  time  level  t*  =  {t-\-  At)  which  can  not  be  computed  using  (  3.10).  The  seaward 
boundary  algorithm  needs  to  be  developed  for  the  cases  of  IWAVE=  1  and  2  where 
the  incident  wave  profile  77,(t)  at  a:  =  0  is  specified  as  input  and  for  IWAVE=3 
where  the  total  free  surface  profile  rj{t)  at  a;  =  0  is  specified  as  input. 
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In  order  to  derive  approximate  seaward  boundary  conditions  for  h  and  q, 
(2.18)  and  (2.20)  are  expressed  in  the  following  characteristic  forms: 


da 

dt 


da  ,  da 


Tb  1  dm 
h  h  dx 


dx 

along  -j-  =  U  +  c 
dt 


dt 


dt 


+  iu 


h  dx 


along 


dx 

dt 


U-c 


(3.24) 

(3.25) 


with 

c  =  y/h  ;  a  —  U  +  2c  ;  ^  =  —U  +  2c  (3.26) 

in  which  c  is  the  normalized  phase  velocity,  whereas  a  and  ^  are  the  characteristic 
variables. 

Assuming  that  U  <  c  and  the  flow  is  subcritical  in  the  vicinity  of  the  seaward 
boundary  where  the  normalized  water  depth  below  SWL  is  dt,  a  and  jd  represent 
the  characteristics  advancing  landward  and  seaward,  respectively,  in  the  vicinity  of 
the  seaward  boundary.  The  total  water  depth  at  the  seaward  boundary  is  expressed 
in  the  form  (Kobayashi  et  al.  1987). 


h{t)  =  dt  +  r){t)  at  X  =  0 


(3.27) 


with 

rf{t)  —  r]i{t)  +  r]r{t)  at  x  =  0  (3.28) 

where  rji  and  rjr  are  the  free  surface  elevations  normalized  by  if'  at  x  =  0  due  to  the 
incident  and  reflected  waves,  respectively.  The  incident  wave  train  may  be  specified 
by  prescribing  the  variation  of  rji  with  respect  to  t  iov  t  >  0.  Alternatively,  the  free 
surface  elevation  measured  at  x  =  0  may  be  specified  as  input  by  prescribing  the 
variation  of  rj  with  respect  to  t  for  t  >  0.  The  normalized  reflected  wave  train 
is  approximately  expressed  in  terms  of  the  seaward  advancing  characteristic  ^  at 

X  =  0 


r)r{t)  ~  -dt-Ct  at  X  =  0 


(3.29) 


where  ^  is  obtained  using  (  3.25)  and  linear  long  wave  theory  is  used  to  derive 
(  3.29).  For  h{t)  at  a:  =  0  calculated  using  (  3.27),  U  =  {2'/h  —  using  (  3.26)  and 
q  =  hU. 

The  correction  term  Ct  in  (  3.29)  introduced  by  Kobayashi  et  al.  (1989)  to 
predict  wave  set-down  and  setup  on  a  beach  may  be  expressed  as 

at  x  =  0  (3.30) 

I  a 


For  incident  regular  waves  on  gentle  slopes,  Ct  may  be  estimated  by  (Kobayashi  et 
al.  1989) 

Ct  =  for  gentle  slopes  (3.31) 

IQdt 

where  the  assumptions  of  linear  long  wave  and  negligible  wave  reflection  were  made 
in  (  3.30)  to  derive  (  3.31).  For  coastal  structures,  wave  reflection  may  not  be 
negligible  but  the  location  of  the  seaward  boundary  may  be  chosen  such  that  Ct  —  0 
on  the  basis  of  (  3.31).  For  incident  irregular  waves,  (  3.31)  may  still  be  used  as 
a  first  approximation  to  improve  the  prediction  of  wave  set-down  and  setup  on  a 
beach.  When  the  measured  time  series  of  77 (t)  at  a;  =  0  specified  as  input  includes 
the  wave  set-down  or  setup  at  a;  =  0,  the  reflected  wave  train  Tir{t)  is  computed 
using  (  3.29)  with  Ct  =  0. 

On  the  other  hand,  the  value  of  m  at  the  seaward  boundary  needs  to  be 
found  using  (2.31).  The  initial  condition  for  m  is  specified  as  m  =  0  at  t  =  0  in 
the  computation  domain  a;  >  0.  The  value  of  m  at  x  =  0  might  be  taken  as  m  =  0 
at  X  =  0  if  the  seaward  boundary  is  located  outside  the  surf  zone.  This  is  because 
the  vertical  variation  of  the  horizontal  velocity  assumed  in  (2.33)  is  caused  by  wave 
breaking  in  this  numerical  model  for  shallow  water  waves.  However,  the  boundary 
condition  of  m  =  0  at  x  =  0  will  yield  m  =  0  for  7  >  0  and  x  >  0  because  m  =  0 
is  a  trivial  solution  of  (2.31).  It  is  hence  required  to  introduce  m  >  0  at  x  =  0  so 
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that  m  >  0  for  t  >  0  and  x  >  0.  One  option  is  to  rewrite  (2.31)  in  terms  of  Uj,  using 
(2.37),  (2.38)  and  (2.39) 


diib  d  CsUb  lub  dh  ^dub\  Tb  +  CBe\ub\ub 

dt  dx  ^  ~  2C2  \h  dx^  dx )  C2h 


with 


CBt  =  CbCJc 


(3.32) 


(3.33) 


in  which  u  is  given  by  (2.28)  with  Ub  =  {U  +  Ub)  and  Ufc  =  0  is  not  a  trivial  solution 
of  (  3.32).  The  value  of  m  =  C2hul  at  x  =  0  may  be  obtained  using  the  value  of  Ub 
at  X  =  0  computed  using  (  3.32)  as  explained  in  the  following. 


3.5.1  Seaward  Boundary  Algorithm  for  Specified  Incident  Regular  and 
Irregular  Waves 

An  explicit  first-order  finite  difference  equation  corresponding  to  (  3.25)  is 
used  to  find  the  value  of  at  x  =  0  and  the  next  time  t*  for  the  cases  of  specified 
incident  regular  and  irregular  waves 


K  =  A  -  -  ciXA  -  A)  +  A( 


^1  + 


(T-fc)! 


h 


where  =  {—U\  -f  2ci)  and  /?2  =  {—U2  +  2C2).  The  right  hand  side  of  (  3.34)  can  be 
computed  for  the  known  values  of  Uj  with  j  =  1  and  2  at  the  present  time  t  where 
the  spatial  nodes  are  located  at  x  =  (j  —  l)Ax.  The  value  of  r)*  at  the  time  t*  is 
calculated  using  (  3.29).  The  incident  wave  profile  specified  as  input  together 
with  (  3.27)  and  (  3.28)  yields  the  value  of  while  t/^  =  —  /3J']  using  the 

definition  of  ^  given  in  (  3.26).  Thus,  the  values  of  h\,  t/f ,  and  =  U^h\  at  x  =  0 
and  the  time  t*  are  obtained. 


As  for  the  value  of  m*  at  x  =  0  and  the  next  time  t*,  an  explicit  first-order 
finite  difference  approximation  of  (  3.32)  is  used  to  obtain  the  value  of  as 
follows: 


(5,): = ^  |£/2  ink  -  u,  (it).]  - 


C3{^b)l 

2Ax 


4-3  {ub)2]  +  K  ^  [(^6)i  +  Cbi  I  iub)i  I  («6)i]}  (3.35) 


The  value  of  is  then  calculated  using  (2.37) 


C2K 


(3.36) 


3.5.2  Seaward  Boundary  Algorithm  for  Specification  of  Total  Incident 
and  Reflected  Waves 


When  the  total  free  surface  comprised  of  both  incident  and  reflected  waves 
7/(t)  at  X  =  0  in  (  3.27)  is  specified  as  input,  the  values  of  rfl  and  at  x  =  0  and 
the  time  t*  are  known.  The  value  oi  at  x  =  0  and  the  time  t*  is  computed  using 
(  3.25)  for  the  characteristic  variable  (3  advancing  seaward  from  the  computation 
domain. 

A  simple  first-order  finite  difference  approximation  of  ( 3.25)  along  the  straight 
line,  dxjdt  =  [U^  —  c\)  <  0,  originating  from  the  point  at  node  1  and  the  time 
t*  =  (t  -{■  At)  may  be  expressed  as 


/^r  —  ^12 + At 


^1  + 


{n)i 


hi 


At  m2  —  mi 
Ax  hi 


(3.37) 


where  /?i2  is  the  value  of  ^  at  the  time  t  and  at  the  location  of  x  =  6x  given  by 


Sx=-{Ui-  ct)  At  >0  (3.38) 

The  numerical  stability  criterion  given  by  (  3.11)  requires  that  Sx  <  Ax. 
As  a  result,  the  point  of  x  =  ^x  is  located  between  nodes  1  and  2.  The  linear 
interpolation  between  the  known  values  of  j3i  and  ^2  at  the  time  t  yields 

Sx 

P12  =  ^1  +  {^2  —  /5i) 


(3.39) 


Using  (  3.38)  and  (  3.39),  (  3.37)  may  be  rewritten  as 


K  =  A  -  {u;  -  (A  -  A)  +  At 


A  + 


(^6)i 


which  corresponds  to  (  3.34)  except  that  [Ui  —  Ci)  in  (  3.34)  is  replaced  by  (C/j  —  c^). 
Eq.  (  3.40)  for  =  (2c*  —Ui)  is  an  implicit  scheme  for  Ui  for  the  known  value  of 
Cj  =  ^(h*).  Solving  (  3.40)  for  yields 


f/i*=  1 


At 

^(A-A) 


{2c: -A 


At 
Ax  L 


c*AP2-I3i)  + 


At 


m2  —  mi 


hi 


0i  + 


in)i 


hi 


(3.41) 


If  the  absolute  value  of  the  denominator  on  the  right  hand  side  of  (  3.41)  becomes 
almost  zero,  this  implicit  algorithm  may  not  be  appropriate.  This  problem  has  never 
been  encountered  so  far  partly  because  the  numerical  stability  criterion  expressed 
as  (  3.11)  generally  requires  a  value  of  At/ Ax  that  is  much  less  than  unity. 

After  Ui  is  computed  using  (  3.41),  the  value  of  is  obtained  from  ^i  = 
(2ci  —  Ui).  The  value  of  ij*  for  the  reflected  wave  profile  is  then  calculated  using 
(  3.29)  at  the  time  t*.  The  value  of  rf*  for  the  incident  wave  profile  is  obtained  from 
Vi  —  {Vi  ~  Vr)  based  on  (  3.28).  The  value  of  m*  is  computed  using  (  3.36)  with 
(  3.35). 


3.5.3  Wave  Reflection  Coefficient 

The  average  reflection  coefficient  r  for  regular  and  irregular  waves  may  be 
estimated  using  the  root-mean-square  values  of  the  time  series  of  rjr  and  rji  as  defined 
by  (  3.23) 

^^iVr)mts/iVi)rms  (3-42) 

which  is  equal  to  the  square  root  of  the  ratio  between  the  time-averaged  reflected 
wave  energy  as  compared  to  the  time-averaged  incident  wave  energy  on  the  basis 
of  linear  wave  theory.  The  reflection  coefficient  as  a  function  of  the  frequency  for 
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irregular  waves  can  be  calculated  using  the  reflected  and  incident  wave  spectra 
computed  from  the  time  series  of  rjr  and  rn  {e.g.,  Kobayashi  et  al.  1990). 


3.6  Wave  Runup 

For  the  case  of  no  wave  overtopping,  the  landward  boundary  of  the  numer¬ 
ical  model  is  located  at  the  moving  shoreline  on  the  slope  where  the  water  depth 
is  essentially  zero.  The  kinematic  boundary  condition  requires  that  the  horizontal 
shoreline  velocity  be  the  same  as  the  horizontal  fluid  velocity.  In  reality,  it  is  difficult 
to  pinpoint  the  exact  location  of  the  moving  shoreline  on  the  slope.  For  the  compu¬ 
tation,  the  shoreline  is  defined  as  the  location  where  the  normalized  instantaneous 
water  depth  equals  a  small  value  8  such  as  ^  =  10“^  as  explained  in  Section  3.1. 

The  following  numerical  procedure  dealing  with  the  moving  shoreline  located 
at  ft  =  ^  is  used  to  obtain  the  values  of  Uj  at  the  next  time  t*  =  {t  At)  for  the 
nodes  j  >  s  where  s  =  integer  indicating  the  wet  node  next  to  the  moving  shoreline 
at  the  present  time  t  such  that  ft^+i  <  8  <  hg.  It  is  noted  that  the  procedure  is 
somewhat  intuitive  and  may  be  improved  since  the  moving  shoreline  tends  to  cause 
numerical  instability. 

1.  After  computing  with  i  =  2,  3,  . . . ,  s  using  (  3.10),  it  is  checked  whether 
ft*_i  <  8,  which  may  be  encountered  during  a  downrush.  This  is  considered  a 
computation  failure  since  the  shoreline  should  not  move  more  than  Ax  because 
of  the  numerical  stability  criterion  of  the  adopted  explicit  method  given  by 
(  3.U). 

2.  It  />;  >  use  a;  =  -  A;_2),  and  u;  =  -  f/.lj),  so  that  the 

water  depth  near  the  shoreline  decreases  landward.  The  following  adjustments 
are  made 

.  if\u:\>\u:_,\,setu:  =  o.9u:_,-, 
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•  ii  h*  <  0,  set  h*  =  0.5^*_i; 

•  and  if  h*  >  set  h*  =  0.9h*_-^. 

Then,  obtain  q*  =  h*U*  based  on  the  adjusted  values  of  h*  and  U*. 

3.  If  h*  <  6,  set  s*  =  (s  —  1)  and  go  to  Step  11.  The  integer  s*  indicates  the  wet 
node  next  to  the  shoreline  at  the  next  time  t*. 

4.  It  h"  >  S,  compute  =  (2h",  -  U;^.,  =  {2U;  -  and  gj,.,  = 

5.  If  <  d,  set  s*  =  s  and  go  to  Step  11. 

6.  If  >  6,  compute  U**  at  the  time  t**  =  {t*  +  At)  using  (  3.10)  with  m  =  0 

where  and  in  (  3.10)  are  replaced  by  U**  and  U*,  respectively.  The 
vertical  variation  of  the  horizontal  velocity  may  be  assumed  to  be  small  in  the 
vicinity  of  the  moving  shoreline.  Improve  the  linearly  extrapolated  values  in 
Step  4  using  the  following  finite  difference  equations  derived  from  (2.18)  and 
(2.20)  with  m  =  0  at  Tfc  =  0: 

A  7* 

?:+!  =  Cl  -  ^  (C  -  hs)  (3.43) 

1  r 

u;+,  =  UU  -  ^  (C."  -  u.)  +  -  hU  +  2Ax  S.j  (3.44) 

The  upper  limit  of  the  absolute  value  of  in  (  3.44)  is  taken  as  8~^  to 

avoid  the  division  by  the  very  small  value.  Calculate  h*^-y  =  qs+i/U*^^. 

7.  If  \U*^i  \  <  6,  set  s*  =  s  and  go  to  Step  11. 

8.  If  <  h*  and  <  8,  set  s*  =  s  and  go  to  Step  11. 

9.  If  <  h*  and  >  8,  set  s*  =  (s  +  1)  and  go  to  Step  11. 
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10.  If  >  h*,  the  linearly  extrapolated  values  of  and  5*^.2  in  Step  4 

are  adopted  in  the  following  instead  of  those  computed  in  Step  6.  Furthermore, 
set  s*  =  (s  + 1)  if  <  hi  and  U*^i  >  8  and  set  s*  =  s  otherwise  where 
and  are  the  adopted  values. 

11.  After  $*  is  obtained,  set  hj  =  0,  Uj  =  0,  gj  =  0  and  rrij  =  0  for  j  >  (s*  +  1) 
since  no  water  is  present  above  the  computational  shoreline.  If  s*  =  (s  +  1), 
set  =  0.  It  is  noted  that  the  smoothing  procedure  given  by  (  3.14)  does 
not  affect  this  shoreline  algorithm. 

Once  the  normalized  water  depth  h  at  the  given  time  is  known  as  a  function  of 
X,  the  normalized  free  surface  elevation,  Zr  =  Z'^fH',  where  the  physical  water  depth 
equals  a  specified  value  can  be  computed  as  long  as  8r  =  {SH H')  >  8.  The  use  of 
the  physical  depth  8[  is  related  to  the  use  of  a  runup  wire  to  measure  the  shoreline 
oscillation  on  the  slope  (e.^'.,  Raubenheimer  et  al.  1995).  The  specified  depth  8^  can 
be  regarded  as  the  vertical  distance  between  the  runup  wire  and  the  slope,  while  the 
corresponding  elevation  Z'  is  the  elevation  above  SWL  of  the  intersection  between 
the  runup  wire  and  the  free  surface.  The  computed  oscillations  of  Zr(t)  for  different 
values  of  8^  can  be  used  to  examine  the  sensitivity  to  8^  of  wave  runup  and  run¬ 
down,  which  are  normally  defined  as  the  maximum  and  minimum  elevations  relative 
to  SWL  reached  by  uprushing  and  downrushing  water  on  the  slope,  respectively.  The 
normalized  runup  R,  run-down  and  setup  Zr  as  well  as  the  root-mean-square 
value  (standard  deviation)  of  Zr{t)  for  given  8^  are  obtained  from  the  computed 
oscillation  of  Zr{t). 


Chapter  4 


COMPARISON  OF  MODEL  VBREAK  TO  DATA  AND 

PREVIOUS  MODEL 

4.1  Comparison  of  Model  to  Previously  Developed  One  Dimensional 
Model 

To  demonstrate  the  effectiveness  of  the  MacCormack  method  in  the  numer¬ 
ical  solutions  of  the  continuity  and  momentum  equations,  the  model  VBREAK  is 
reduced  to  a  one-dimensional  model  and  compared  with  a  previously  developed  one¬ 
dimensional  model.  VBREAK  implements  the  explicit  first  order  finite  difference 
of  (  3.32)  to  obtain  the  value  of  and  thus  m*,  the  momentum  correction  factor 
at  the  next  time  level  at  the  seaward  boundary.  When  the  value  of  m\  is  set  to  zero 
at  all  times,  the  computed  values  of  m*j  at  any  j  are  zero  for  all  times  everywhere 
inside  the  computational  domain.  Thus  the  computed  horizontal  velocity  has  no 
vertical  variation.  For  the  entirety  of  this  Section  4.1,  VBREAK  is  reduced  to  a 
one-dimensional  model  through  the  specification  of  zero  vertical  variation  at  the 
seaward  boundary. 

The  model  VBREAK  as  a  a  one-dimensional  model  is  compared  with  the 
previously  developed  IBREAK;  the  two  models  employ  different  numerical  methods 
to  solve  the  same  governing  equations.  Use  is  made  of  the  1:2.5  riprap  revetment 
test  conducted  by  Ahrens  (1975).  The  comparison  of  VBREAK  to  IBREAK  is  the 
only  example  of  a  steep  slope  computed  in  this  report.  The  seaward  boundary  for 
the  computation  is  taken  at  the  toe  of  the  1:2.5  slope  where  the  water  depth  below 


SWL  d[  is  15  ft.  Cnoidal  waves  of  8.5  s  period  with  a  height  of  3.06  ft  are  the  time 
series  of  the  incident  wave  at  the  seaward  boundary.  A  friction  factor  of  0.3  for  the 
riprap  slope  (Kobayashi  et  al.  1987)  and  a  Courant  number  of  0.4  are  used  for  the 
computation.  Although  the  numerical  model  VBREAK  is  found  to  be  stable  for 
Courant  numbers  up  to  0.7,  the  value  Cn  =  0.4  is  used  so  as  to  be  consistent  with 
other  comparisons  throughout  this  report.  A  numerical  damping  coefficient  k  of  0.1 
is  sufficient  for  damping  high  frequency  oscillations  throughout  the  six  wave  periods 
from  t  =  0  to  t  =  tmax  =  6.  The  steep  slope  case  suffers  less  numerical  oscillation 
problem  at  the  rear  of  the  steep  front  than  the  gently  sloped  cases  computed  in 
the  subsequent  sections.  Tables  4.1  and  4.2  are  the  primary  input  and  output 
from  the  model  VBREAK  as  explained  by  Kobayashi  and  Johnson  (1995).  The 
resultant  dimensionless  parameters  for  this  case  are;  surf  similarity  parameter  ^  = 
4.4;  the  ratio  of  horizontal  to  vertical  length  scales  a  =  27.6;  Ursell  number  =  30.1. 
As  evident  in  the  primary  output  file,  the  reflection  coefficient  (r  =  0.55)  is  large 
as  expected  for  the  steep  slope  and  surf  similarity  parameter  ^  =  4.4.  Although 
information  on  the  cubic  profile  parameter  a  is  included,  it  is  not  utilized  in  the 
computations  until  the  following  sections. 

Figure  4.1  depicts  the  six  Cnoidal  waves  used  as  input  to  the  models  VBREAK 
and  IBREAK  as  well  as  the  reflected  wave.  Incident  and  reflected  waves  are  found 
to  be  essentially  the  same  for  both  models.  The  reflected  wave  is  of  considerable  size 
relative  to  the  incident  wave  and  periodicity  is  reached  within  two  wave  periods. 

The  runup  for  the  computational  shoreline  based  on  the  water  depth  = 
0.787  in  is  shown  in  Figure  4.2  where  the  incident  wave  train  arrives  at  the  shoreline 
at  t  ~  0.5.  The  computed  runup  Z  is  the  normalized  shoreline  elevation  above  the 
SWL  and  is  plotted  against  the  normalized  time.  Following  an  initial  depression 
of  the  free  surface,  the  runup  consists  of  the  setup  ^  ~  1.0  above  the  SWL  and  a 
periodic  component  of  amplitude  ~  1.0.  Figure  4.2  shows  that  the  two  models  yield 


37 


Table  4.1:  Primary  input  data  file  for  VBREAK  in  comparison  with  IBREAK  . 


=  number  of  comment  lines 


Ahren  (1975)  Test  12 


<— IWAVE 

<— IBOT 

<~INCLCT 

<“IENERG 

<~ITEMVA 

<“ISPAEU 


5. 

6. 

<— TSTAT,TMAX 

100 

<— STILL 

0.001 

0.400 

0.1  <--DELTA 

3.06 

8.5 

<— HREF,TREF 

1.0 

<— KS 

24001 

<— NPINP 

4.0 

.10 

<-“APROFL,CMIXL 

15. 

0.4 

<— DSEAP,SLSURF 

1 

<— NBSEG 

60.000 

9.000 

.30 

1201 

<— NPOUT 

3 

0 . 157480 

0.787402 

1.574803 

<— NDELR 

4 

<-’-N0N0DS 

91  101 

111  121 

<— NODLOC(l)  (2) 

<— NOTIML 
5.75 


<— TIMSPA 
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Table  4.2:  Primary  output  file  for  VBREAK  in  comparison  with  IBREAK. 


Ahren  (1975)  Test  12 


WAVE  CONDITION 


Cnoidal  Incident  Wave  at  Seaward  Boundary 
Normalized  wave  height  KS  =  1.000000 

1-m  =  0.112815261D+00 

E  =  0.111509461D+01 

K  =  0.252149782D+01 


Reference  Wave  Period 


8.500000  sec. 


Reference  Wave  Height 
Depth  at  Seaward  Boundary 
Norm.  Depth  at  Seaw.  Bdr. 


3.060000  feet 
15.000000  feet 
4.902 


Included  Correction  Term  CT 
0  =  no;  1  =  yes  INCLCT  =  0 

Normalized  Wave  Length  =  12.144 

"Sigma"  =  27.573 

Ursell  Number  =  30.084 

Surf  Similarity  Parameter  =  4.400 

Input  Wave  Train  from  Time=0  to  TMAX 
Computed  or  Read  at  Normalized  Rate  DELTI  = 


0.000250 


Parameters  of  Vertical  Velocity  Vaoriations 


Cubic  Profile  Parameter  APROFL  = 
Mixing  Length  Parameter  CMIXL  = 
Momentum  Flux  Coefficient  C2  = 
Kinetic  Energy  Flux  Coeff.  C3  = 
Energy  Dissipation  Coeff.  CB  = 
Coefficient  of  DB  CBL  = 


4.000000 

0.100000 

0.485714 

0.000000 

12.342857 

3.403313 


BOTTOM  GEOMETRY 

Norm,  Horiz.  Length  of 

Computation  Domain  =  0.711121 

Number  of  Segments  =  1 


SEGMENT  XBSEG(I)  ZBSEG(I)  BFFSEG(I) 

I  feet  feet 


.X=0 

1 


0.000000  “15.000000 

60.000000  9.000000 


0.300000 
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Table  4.2:  -  Continued 


COMPUTATION  PARAMETERS 


Normalized  DX  = 

Normalized  DELTA  = 

C our ant  Number  = 

Must  not  exceed  unity 
Numerical  Damping  Coefficient  = 

Must  be  zero  or  positive 
Normalized  Computation  Duration  TMAX  = 
Statistical  Calculations  Start 
when  Time  is  equal  to  TSTAT 

Total  Number  of  Spatial  Nodes  JMAX  = 
Number  of  Nodes  Along  Bottom  Below  SWL 

STILL  = 


0. 44445  lD-‘02 
O.lOOOOOE-02 
0.400 

0.1000 

6.000000 

5.000000 

161 

100 


Storing  Temporal  Variations  from  Time  =  0 

to  TMAX  at  Normalized  Rate  DELTO  = 

Wave  Runup  Time  Series  Stored  for 

NDER  = 

Time  Series  of  ETA,  U,  and  UB 

Stored  at  NONODS  = 

Spacial  Variations  of  ETA,  U,  and  UB 

Stored  at  NOTIML  = 


0.005000 
3  Water  Depths 
4  Nodes 
5  Time  Levels 


Maximum  time  step  =  0.82969E-03 

Minimiun  time  step  =  0.50928E-03 


REFLECTION  COEFFICIENT 


ETARRMS/ETAIRMS  =  0.550 

INCIDENT  AND  REFLECTED  WAVES 


Max 

Min 

Mean 

RMS 

Inc . 

0.6287 

-0.3713 

0.0000 

0.3472 

Ref. 

0.2790 

-0.2720 

0.0143 

0.1910 

SHORELINE 

OSCILLATIONS 

Largest  Node  Number  Reached  by  Computational  Shoreline 

SMAX 

=  138 

I 

DELTAR(I) 

RUNUP (I) 

RUNDOWN (I) 

SETUP(I)  RMS(I) 

[inch] 

Ru 

Rd 

Zr  Rrms 

1 

0.157 

1.831 

0.641 

1.295  0.365 

2 

0.787 

1.811 

-0.322 

0.874  0.668 

3 

1.575 

1.803 

-0.746 

0.670  0.813 
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Figure  4.1:  Incident  and  reflected  waves  at  the  seaward  boundary:  Incident ( — ) 
Reflected!-  •  •). 


Runup 


Figure  4.2:  Computed  time  series  of  normalized  runup  with  =  0.787  in  : 
IBREAK(— );  VBREAK(-  -  -). 

results  that  are  almost  indistinguishable.  Two  additional  computational  shoreline 
depths  are  also  specified  but  are  not  shown  here  as  the  results  are  basically  similar. 

Figures  4.3  and  4.4  depict  the  computed  spatial  variation  of  the  normalized 
free  surface  rj  and  the  normalized  depth  averaged  velocity  U  at  five  times  throughout 
the  final  wave  period.  The  plots  for  t  =  5  and  t  =  6  are  identical  because  of  the 
periodicity.  In  each  panel,  the  variations  of  7]  computed  by  VBREAK  and  IBREAK 
match  almost  exactly.  The  planar  bottom  is  shown  in  each  panel  as  a  straight  line. 
On  the  other  hand,  the  spatial  variations  of  the  depth  averaged  horizontal  velocities 
computed  by  the  two  models  differ  only  slightly.  The  greatest  difference  occurs  at 
the  shoreline  where  the  slight  irregularity  of  the  results  of  IBREAK  seems  physically 
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implausible. 

The  cross-shore  variations  of  the  maximum,  minimum,  and  mean  of  the  nor¬ 
malized  free  surface  with  the  planar  bottom  profile  and  the  depth  averaged  velocity 
are  depicted  in  Figure  4.5.  The  free  surface  statistics  are  practically  the  same.  The 
wave  height  increases  as  the  wave  shoals  and  is  reflected  from  the  steep  slope.  The 
mean  free  surface  reveals  a  small  set-down  before  wave  uprush  occurs  at  approxi¬ 
mately  X  =  0.3.  A  significant  setup  occurs  at  the  shoreline  as  shown  in  Table  4.2. 
The  depth  averaged  velocity  statistics  are,  again,  essentially  the  same  for  the  two 
numerical  models.  A  small  discrepancy  is  apparent  in  the  maximum  depth  averaged 
velocity  near  the  shore  as  the  model  VBREAK  predicts  a  value  that  is  approxi¬ 
mately  20%  larger  than  the  prediction  of  IBREAK.  The  mean  of  the  horizontal 
velocity  U  is  a  small  negative  value  throughout  the  domain.  This  is  a  return  current 
caused  by  forward  velocities  under  the  crest  with  larger  water  depth. 

The  cross-shore  variation  of  the  mean  volume  flux  is  cited  as  an  indicator 
of  the  computationah  accuracy  (Kobayashi  et  al.  1989).  Theoretically,  the  no  flux 
condition  into  the  impermeable  beach  would  dictate  a  net  volume  flux  of  zero.  The 
time  averaged  continuity  equation  corresponding  to  (2.18)  indicates  that  the  time 
averaged  volume  flux  would  be  zero  throughout  the  domain.  Figure  4.6  demon¬ 
strates  that  the  flux  is  indeed  very  much  less  than  one  for  both  models.  IBREAK 
demonstrates  a  greater  accuracy  than  the  developed  model  for  this  case.  However, 
the  net  volume  flux  computed  by  VBREAK  of  approximately  0.003  or  less  is  still 
acceptably  small. 

The  normalized  energy  quantities  for  both  models  are  shown  in  Figure  4.7. 
The  normalized  mean  specific  wave  energy  shown  in  the  first  panel  increases  as 
the  waves  shoal  then  falls  to  zero  at  the  shoreline.  The  normalized  mean  energy 
flux  as  defined  in  Chapter  2  falls  monotonically  to  zero  at  the  shoreline  as  energy  is 
dissipated  through  bottom  friction  and  wave  breaking.  Panel  three  shows  the  energy 
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Figure  4.6:  Computed  cross-shore  variation  of  normalized,  time  averaged  volume 
flux:  IBREAK(— );  VBREAK(-  - 


dissipated  in  the  bottom  boundary  layer  as  explained  in  Section  2.  For  the  non¬ 
breaking  waves  on  the  steep  riprap  slope,  the  wave  energy  dissipation  by  the  bottom 
friction  is  dominant.  Panel  four  indicates  that  VBREAK  dissipates  only  a  small 
fraction  of  energy  in  breaking  while  IBREAK  does  not  include  energy  dissipation 
due  to  the  vertical  shear.  This  is  expected  as  the  breaking  dissipation  is  dependent 
on  a  vertical  structure  as  addressed  in  Section  2.2.  The  final  panel  depicts  the 
numerical  dissipation  from  both  models.  The  numerical  dissipation  is  calculated 
through  an  energy  balance  as  outlined  in  (2.46)  and  (2.47)  where  discrepancies  are 
attributed  to  numerical  dissipation  related  to  no  explicit  wave  breaking  modeling 
in  IBREAK  and  the  deficiency  of  the  simplified  wave  breaking  modeling  adopted 
in  VBREAK.  It  should  be  noted  that  the  numerical  dissipation  is  small  relative  to 
the  physical  dissipation  due  to  the  bottom  friction. 

4.2  Comparison  of  Model  to  Data  of  Stive  (1980) 

The  model  VBREAK  is  compared  with  the  comprehensive  measurements  of 
test  1  presented  by  Stive  (1980)  and  Stive  and  Wind  (1982).  Because  the  numerical 
model  VBREAK  predicts  the  vertical  variations  of  the  horizontal  velocity,  the  com¬ 
parison  of  the  measured  and  computed  velocities  can  be  made  without  any  ambigu¬ 
ity.  Likewise,  the  free  surface  characteristics  might  demonstrate  greater  agreement 
with  the  measured  data  as  a  result  of  a  realistic  vertical  structure.  In  Stive’s  test 

l,  the  incident  regular  waves  with  the  period  T'  =  1.79  s  broke  as  spilling  breakers 
on  the  1:40  concrete  beach.  The  seaward  boundary  for  the  computation  is  taken 
to  be  at  the  still  water  depth  d^  =  0.2375  m,  where  the  near-breaking  wave  profile 
was  shown  to  be  similar  to  the  cnoidal  wave  profile  as  explained  by  Kobayashi  et 
al.  (1989).  The  measured  wave  height  at  the  seaward  boundary  was  H'  =  0.172 

m.  The  free  surface  and  velocity  characteristics  have  been  found  to  be  relatively 
insensitive  to  the  choice  of  the  cubic  velocity  profile  parameter,  a.  In  this  particular 
comparison,  the  parameter  a  is  set  to  3.0.  The  friction  factor  of  0.05  is  used  as  in 
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Figure  4.7:  Computed  cross-shore  variation  of  normalized  specific  wave  energy, 
energy  flux,  bottom  dissipation,  breaking  dissipation,  and  numerical 
dissipation:  IBREAK(— );  VBREAKf-  - 


the  previous  computation  by  Kobayashi  et  al.  (1989).  A  Courant  number  of  0.3  is 
adopted  for  the  stable  computation.  Numerical  stability  is  more  difficult  to  maintain 
during  the  longer  computation  for  the  gentle  slope  of  Stive  compared  with  the  steep 
slope  case  in  Section  4.1.  The  numerical  damping  coefficient  k  is  taken  to  be  1.0 
in  an  effort  to  suppress  the  high  frequency  oscillations  that  occur  at  the  rear  of  the 
breaking  wave  crest.  Tables  4.3  and  4.4  represent  the  primary  input  and  output  files 
of  the  numerical  model  VBREAK.  The  resultant  dimensionless  parameters  were: 
surf  similarity  parameter  ^  =  0.135;  the  ratio  of  horizontal  to  vertical  length  scales 
a  —  13.5;  and  Ursell  number  =  121.7. 

The  normalized  incident  wave  specified  as  input  to  the  model  and  the  com¬ 
puted  reflected  wave  are  depicted  in  Figure  4.8.  The  reflection  coefficient  r  =  0.008 
is  very  small  as  expected  from  the  small  surf  similarity  parameter  ^  =  0.135. 

With  a  computational  shoreline  based  on  the  water  depth  =  1  cm,  the 
computed  runup  is  shown  in  Figure  4.9  where  the  incident  wave  train  arrives  at  the 
shoreline  at  f  ~  5.  The  periodic  runup  consists  of  the  steady  setup  fj  ~  0.06  from 
the  SWL  and  a  relatively  small  oscillatory  component  of  amplitude  0.015.  The 
swash  is  much  smaller  than  the  swash  on  the  steep  slope  presented  in  Section  4.1. 

Figures  4.10  and  4.11  depict  the  computed  spatial  variation  of  the  normal¬ 
ized  free  surface  i]  and  the  depth  averaged  velocity  U  at  five  times  throughout  the 
final  wave  period.  The  plots  for  t  =  29  and  t  =  30  are  identical  because  of  the 
periodicity.  As  the  wave  propagates  shoreward,  the  saw  tooth  profile  develops.  The 
high  frequency  numerical  oscillations  are  apparent  following  the  crest  of  the  wave; 
The  numerical  smoothing  described  in  Section  3.2  is  intended  to  reduce  these  oscil¬ 
lations.  The  previous  one-dimensional  models  based  on  the  Lax-Wendroff  scheme 
also  developed  these  unrealistic  numerical  irregularities.  Figure  4.11  also  shows 
the  near  bottom  velocity  calculated  by  model  VBREAK.  The  horizontal  velocity 
Ub  =  (U  -{■  Ub)  immediately  outside  the  bottom  boundary  layer  is  attenuated  due  to 


Table  4.3:  Primary  input  data  file  for  Stive’s  test  1. 


3 


=  number  of  comment  lines 


Stive  1980  Test  1 


1 

1 

1 

1 

1 

1 

1 

29. 

190 

0.001 

.172 

1.0 

9001 

3.0 

.2375 

1 

18.000 

3001 

1 

1. 

6 

1  41 

141 

5 

29.0 


30. 


0.300 

1.79 


.10 

0.025 

.025  .05 


61  81  101 

29.25 


1.0 


<-“ISYST 
<--IWAVE 
<— IBOT 
<~~INCLCT 
<— lENERG 
<— ITEMVA 
<-- ‘ISPAEU 
<— TSTAT,TMAX 
<— STILL 

<— DELTA ,  COURNO  .DKAPPA 

<— HREF,TREF 

<— KS 

<— NPINP 

<— APROFL,CMIXL 

<— DSEAP,SLSURF 

<— NBSEG 

<— WBSEG  ( 1 )  ,  TBSLOP  ( 1 )  ,  BFFSEG  ( 1 ) 

<“-NP0UT 

<-~NDELR 

<— DELRP(l) 

<— NONODS 

<— N0DL0C(1,2,3,4,5) 

<— N0DL0C(6) 

<— NOTIML 
30.00 


29.50 


29.75 


<--TIMSPA(l,2,3,4,5) 


Table  4.4:  Primary  output  file  for  Stive’s  test  1. 


Stive  1980  Test  1 


WAVE  CONDITION 


Cnoidal  Incident  Wave  at  Seaward  Boundary 
Normalized  wave  height  KS  = 

1-m  =  0.113153458D-02 

E  =  0.100242146D+01 

K  =  0.477945412D+01 

Reference  Wave  Period 
Reference  Wave  Height 
Depth  at  Seaward  Boundary 
Nonn.  Depth  at  Seaw.  Bdr, 

Included  Correction  Term  CT 
0  =  no;  1  =  yes  INCLCT 

Normalized  Wave  Length 
"Sigma'* 

Ur sell  Number 
Surf  Similarity  Parameter 
Input  Wave  Train  from  Time=0  to 


1.000000 


1.790000  sec. 
0.172000  meters 
0.237500  meters 
1.381 


^  1 

12.963 

13.518 

121.692 

0.135 

TMAX 


Computed  or  Read  at  Normalized  Rate  DELTI  = 


0.003333 


Paxcutieters  of  Vertical  Velocity  Variations 
Cubic  Profile  Parameter  APROFL  =  3.000000 

Mixing  Length  Parameter  CMIXL  =  0.100000 

Momentum  Flux  Coefficient  C2  =  0.548214 

Kinetic  Energy  Flux  Coeff.  C3  =  -0.069420 

Energy  Dissipation  Coeff.  CB  =  15.163393 

Coefficient  of  DB  CBL  =  2,049839 


BOTTOM  GEOMETRY 


Norm.  Horiz.  Length  of 

Computation  Domain 
Number  of  Segments 


7.741422 

1 


SEGMENT  WBSEG(I)  TBSLOP(I)  BFFSEG(I) 

I  meters 


1  18.000000 


0.025000 


0.050000 


Table  4.4:  -  Continued 


COMPUTATION  PARAMETERS 


Normalized  DX  = 

Normalized  DELTA  = 

Courant  Number  = 

Must  not  exceed  unity 
Numerical  Damping  Coefficient  = 

Must  be  zero  or  positive 
Normalized  Computation  Duration  TMAX  = 
Statistical  Calculations  Start 
when  Time  is  equal  to  TSTAT= 

Total  Number  of  Spatial  Nodes  JMAX  = 
Number  of  Nodes  Along  Bottom  Below  SWL 

STILL  = 

Storing  Temporal  Variations  from  Time  = 
to  TMAX  at  Normalized  Rate  DELTO  = 
Wave  Runup  Time  Series  Stored  for 

NDER  = 

Time  Series  of  ETA,  U,  and  UB 

Stored  at  NONODS  = 

Spacial  Variations  of  ETA,  U,  and  UB 

Stored  at  NOTIML  = 


0.215040D-01 

O.lOOOOOE-02 

0.300 

1.0000 

30.000000 

29.000000 

361 

190 

0.010000 
1  Water  Depths 
6  Nodes 
5  Time  Levels 


Maximum  time  step 
Minimum  time  step 


0.73200E-02 

0.27400E-02 


REFLECTION  COEFFICIENT 
ETARRMS/ETAIRMS  =  0.008 

INCIDENT  AND  REFLECTED  WAVES 


Inc . 
Ref. 


Max 

0.7906 

-0.0354 


Min 

-0.2088 
-0 . 0427 


Mean 

0.0000 

-0.0388 


RMS 

0.3096 

0.0024 


SHORELINE  OSCILLATIONS 

Largest  Node  Number  Reached  by  Computational  Shoreline 

SMAX  =  204 


I  DELTAR(I) 
[cm] 


RUNUP (I)  RUNDOWN (I)  SETUP(I)  RMS (I) 
Ru  Rd  Zr  Rrms 


1.000 


0.074 


0.047 


0.059 


0.009 
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Inddent  and  Reflected  Waves  at  Boundary 


Figure  4.8:  Normalized  incident  and  reflected  waves  at  the  seaward  boundary. 
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Figure  4.9:  Computed  normalized  time  series  of  runup  with  —  1  cm. 

bottom  friction.  The  magnitude  of  the  computed  bottom  velocity  is  slightly  smaller 
than  that  of  the  depth  averaged  velocity  because  Ub  is  assumed  to  be  out  of  phase 
with  the  depth  averaged  velocity  U  in  (2.44).  Near  the  bottom,  the  horizontal  ve¬ 
locity  develops  a  kink  as  the  depth  averaged  velocity  U  changes  sign  and  the  value 
of  Ub  changes  abruptly.  This  unrealistic  kink,  which  is  apparent  in  Figure  4.11,  sug¬ 
gests  the  shortcoming  of  the  vertical  variation  of  the  horizontal  velocity  assumed  in 
(2.33). 

Figure  4.12  shows  the  computed  and  measured  maximum  free  surface  eleva¬ 
tions  throughout  the  domain  for  the  final  wave  period.  The  measured  and  computed 
minimum  free  surface  elevations  are  depicted  in  Figure  4.13.  The  wave  height  com¬ 
parison  is  shown  in  Figure  4.14,  defined  cis  the  difference  between  the  maximum  and 
minimum  elevations.  The  data  of  Stive  presented  in  these  plots  is  read  from  Figure 
4  of  Stive  and  Wind  (1982).  The  wave  height  prediction,  appears  to  be  better  in 


55 


Figure  4.11:  Computed  cross-shore  variation  of  normalized  depth  averaged  veloc¬ 
ity  U  and  near  bottom  velocity  uj  at  5  time  levels,  t  =  29.00,  29.25, 
29.50,  29.75,  and  30.00:  U  (-  -  -);  (•  •  •)• 
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Measured  and  Computed  Maximum  eta 


Figure  4.12:  Cross-shore  variation  of  the  measured  and  computed  maximum  nor¬ 
malized  free  surface. 

part,  due  to  the  cancellation  of  the  overprediction  and  underprediction  in  the  com¬ 
puted  maximum  and  minimum  elevations.  The  numerical  model  underpredicts  the 
rapid  decrease  of  the  wave  height  immediately  landward  of  the  break  point  where 
the  wave  height  is  the  largest. 

The  computed  and  measured  time-averaged  free  surface  elevation,  set-down 
or  setup,  is  depicted  in  Figure  4.15.  The  computed  mean  free  surface  is  notably 
different  in  shape  from  that  computed  with  IBREAK  in  Figure  11  of  Kobayashi  et 
al  (1989). 

The  maximum  and  minimum  horizontal  velocities  at  four  spatial  positions 
are  depicted  in  Figure  4.16.  In  each  panel  the  dashed  line  represents  the  computed 
horizontal  velocity  distribution  and  the  solid  line  is  the  computed  depth  averaged 
velocity.  The  vertical  axis  is  the  ratio  of  z I d  where  z  is  the  normalized  vertical 
coordinate  and  d  is  the  normalized  still  water  depth.  This  is  slightly  different  from 
the  vertical  axis  of  Figure  7  in  Stive  (1980)  where  the  vertical  coordinate  equaled 
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Figure  4.13:  Cross-shore  variation  of  the  measured  and  computed  minimum  nor¬ 
malized  free  surface. 
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Figure  4.14:  Cross-shore  variation  of  the  measured  and  computed  normalized  wave 
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Figure  4.15:  Cross-shore  variation  of  the  measured  and  computed  mean  normal¬ 
ized  free  surface. 

zero  at  the  mean  water  level.  It  is  worth  noting  that  the  maximum  and  minimum 
values  in  Figure  4.16  are  obtained  at  each  elevation  ^  without  regard  to  the  vertical 
phase  differences.  Panel  one  indicates  that  although  the  specified  and  measured  free 
surface  elevations  match  at  the  seaward  boundary  reasonably  well,  the  maximum 
horizontal  velocity  is  considerably  overpredicted.  Likewise,  at  the  subsequent  loca¬ 
tions  X  =  1.29,  2.15,  and  3.011  the  predicted  maximum  velocity  is  larger  than  the 
measured.  The  greatest  computed  variation  of  the  velocity  over  the  depth  occurs 
after  breaking,  at  x  =  1.29.  This,  however,  does  not  correspond  well  to  the  mea¬ 
sured  data  below  the  trough  level  that  display  virtually  no  variation  with  depth  at 
X  =  1.29. 

The  model  VBREAK  exhibits  the  small  seaward  oriented  return  current  as 
explained  by  Kobayashi  et  al.  (1989).  Figure  4.17  shows  the  computed  return 
current  and  the  measured  normalized  undertow  below  the  trough  level  at  several 
locations  throughout  the  water  depth.  The  return  current  is  calculated  as  the  time 


Horizontal  Velocity  Data 


average  of  the  horizontal  velocity  at  a  given  elevation.  The  measured  values  are  taken 
from  Figure  8  based  on  Stive’s  data  from  Svendsen  (1984).  The  numerical  model 
predicts  the  order  of  magnitude  of  the  undertow  but  does  not  predict  its  vertical 
variation  well,  perhaps  because  the  model  does  not  account  for  the  landward  mass 
flux  caused  by  a  roller  above  the  trough. 

For  completeness,  the  cross-shore  variations  of  the  maximum,  minimum  and 
mean  of  the  normalized  free  surface,  the  depth  averaged  velocity,  and  the  near 
bottom  velocity  are  depicted  in  Figure  4.18.  The  near  bottom  velocity  is  slightly 
smaller  than  the  depth  averaged  velocity  U  as  expected. 

As  in  Section  4.1,  the  cross-shore  variation  of  the  time- averaged  volume  flux 
is  used  as  an  indicator  of  the  computational  accuracy.  Figure  4.19  demonstrates 
that  the  computed  flux  is  nearly  zero  to  satisfy  the  no  flux  boundary  condition  into 
the  impermeable  beach. 

The  normalized  energy  quantities  are  shown  in  Figure  4.20.  Notable  is  the 
fact  that  the  numerical  dissipation  dominates  over  the  physically  tenable  mecha¬ 
nisms  such  as  breaking  dissipation  and  bottom  dissipation  included  explicitly  in 
VBREAK  .  The  value  of  the  cubic  velocity  profile  parameter  a  has  been  set  to  3.0 
in  order  to  increase  the  breaking  dissipation  as  explained  in  Section  2.2.  However, 
the  breaking  dissipation  is  still  small.  The  increase  of  the  mixing  length  parameter 
Cl  in  2.39  from  0.1  to  0.2  causes  the  decrease  of  resulting  in  little  change  in 
the  breaking  dissipation.  This  indicates  that  the  velocity  profile  assumed  in  (2.33) 
does  not  describe  well  the  wave  energy  dissipation  due  to  breaking  on  the  gentle 
slope.  This  probably  derives  from  the  fact  that  the  wave  front  (roller)  is  not  mod¬ 
eled  specifically. 

The  measured  and  computed  temporal  variations  of  the  free  surface  are  com¬ 
pared  at  four  spatial  locations  throughout  the  domain  as  shown  in  Figure  4.21.  The 
variation  of  the  free  surface  from  the  mean  water  level  fj  for  the  last  wave  from 
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Figure  4.18:  Computed  cross-shore  variation  of  maximum,  minimum,  and  mean 
of  the  normalized  free  surface,  depth  averaged  velocity,  and  near 
bottom  velocity. 
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Figure  4.20:  Computed  cross-shore  variation  of  normalized  specific  wave  energy, 
energy  flux,  bottom  dissipation,  breaking  dissipation,  and  numerical 
dissipation. 


t  =  29.0  to  t  =  30.0  is  shown  in  each  panel.  The  crest  of  the  computed  wave  form 
has  been  matched  by  hand  with  the  measured  crest;  therefore,  Figure  4.21  should 
indicate  only  the  comparison  of  the  predicted  and  measured  wave  shapes.  In  fact, 
differences  between  the  computed  and  actual  phase  speeds  are  expected,  based  on 
the  comparison  of  VBREAK  to  the  data  set  of  Cox  (1995)  which  preserved  the 
phase  information  throughout  the  domain.  Panel  one  shows  the  degree  to  which  the 
specified  incident  cnoidal  wave  agrees  with  the  measured  data  of  Stive.  Again,  the 
high  frequency  oscillations  following  the  wave  crest  are  apparent  in  Figure  4.21. 

4.3  Comparison  of  Model  to  Data  of  Cox  et  al  (1995) 

The  model  VBREAK  is  finally  compared  with  the  data  collected  by  Cox  et  al. 
(1995)  that  included  detailed  velocity  profiles  inside  the  surf  zone.  The  experiment 
was  performed  in  a  wave  flume.  The  1:35  beach  was  constructed  of  Plexiglas  and 
had  a  thin  layer  of  sand  glued  to  the  surface  to  increase  the  bottom  roughness. 
Figure  4.22  depicts  the  experimental  setup  and  measurement  locations.  At  each 
measuring  line,  water  velocity  measurements  were  taken  at  approximately  twenty 
locations  throughout  the  water  column. 

The  constant  wave  period  throughout  the  experiment  was  2.2  s.  The  measur¬ 
ing  line,  L2  was  located  at  the  break  point  defined  by  Cox  et  al  (1995)  as  the  point 
where  aeration  occurred  in  the  tip  of  the  wave  crest.  The  spilling  wave  developed 
into  a  turbulent  bore  landward  of  measuring  line  3. 

The  free  surface  and  velocities  were  measured  for  the  duration  of  the  last 
50  waves  out  of  300  at  each  location.  The  numerical  model  VBREAK  is  run  with 
300  waves  in  order  to  be  consistent  with  the  experimental  procedure.  The  statistical 
measurements  are  based  on  the  last  50  waves  in  the  same  way  as  the  measured  data. 

The  free  surface  and  velocity  characteristics  have  been  found  to  be  insensi¬ 
tive  to  the  choice  of  the  cubic  velocity  profile  parameter  a.  In  the  following,  the 
parameter  a  is  set  to  3.0  as  in  Sections  4.1  and  4.2. 


Measuring  Lines 


Wave  Maker 


Figure  4.22:  Experimental  setup  (Cox  et  al.  1995). 

The  friction  factor  =  0.05  is  the  same  as  that  for  the  gentle  slope  specified 
in  Section  4.2.  An  argument  could  be  made  that  the  friction  factor  for  the  quasi- 
two-dimensional  model  VBREAK  would  need  to  be  recalibrated  because  the  bottom 
friction  is  based  on  a  near  bottom  velocity  Ub  rather  than  a  depth  averaged  velocity 
U.  However,  the  difference  between  the  near  bottom  and  depth  averaged  velocities 
is  not  large  enough  to  warrant  this  effort.  This  is  especially  true  considering  that 
the  bottom  dissipation  is  secondary  in  the  surf  zone  and  that  the  constant  friction 
factor  is  itself  an  approximation. 

A  Courant  number  of  0.4  maintained  numerical  stability  and  has  been  used 
in  both  of  the  following  comparisons. 

4.3.1  Model  Initiation  before  Breaking 

The  model  VBREAK  is  initially  compared  to  the  data  with  the  seaward 
boundary  at  measuring  line  1  depicted  in  Figure  4.22.  The  height  of  the  incident 
regular  waves  at  measuring  line  1  was  H'  —  13.22  cm.  The  corresponding  still  water 
depth  was  28.00  cm.  Tables  4.6  and  4.7  represent  the  primary  input  and  output  from 
the  model  VBREAK.  The  resultant  dimensionless  parameters  were:  surf  similarity 


Table  4.5:  Horizontal  locations  and  still  water  depths  at  six  measuring  lines. 


Line  No. 

(cm) 

d 

(cm) 

LI 

0 

28.00 

L2 

240 

21.14 

L3 

360 

17.71 

L4 

480 

14.29 

L5 

600 

10.86 

L6 

720 

7.43 

parameter  ^  =  0.216;  the  ratio  of  horizontal  to  vertical  length  scales  u  —  19.0;  Ursell 
number  =  74.0. 

The  normalized  incident  wave  profile  measured  at  line  1  is  specified  as  input 
to  VBREAK  and  is  depicted  in  Figure  4.23.  Only  the  last  fifty  waves  are  utilized 
in  the  model  to  data  comparison  and  thus  only  those  fifty  waves  are  depicted  in 
Figure  4.23.  The  total  of  the  300  waves  is  constructed  through  the  repetition  of 
six  sets  of  the  50  waves  depicted  in  Figure  4.23.  The  calculated  reflected  wave  is 
insignificant  relative  to  the  incident  wave.  The  computed  wave  reflection  coefficient 
r  =  0.01  is  consistent  with  the  small  surf  similarity  parameter  ^  =  0.216  for  the  1:35 
slope. 

The  runup  computed  based  on  a  computational  shoreline  of  water  depth 
=  0.5  cm,  is  depicted  in  Figure  4.24.  As  in  the  previous  figure,  only  the  last  50  of 
the  300  waves  are  shown.  The  computed  runup  Z  is  the  shoreline  elevation  above 
the  SWL  and  is  plotted  against  the  normalized  time.  The  incident  wave  shown 
in  Figure  4.23  is  not  breaking  but  is  not  perfectly  periodic.  The  runup  on  the 
gentle  slope  amplifies  the  nonperiodic  components  of  the  incident  waves  shown  in 
Figure  4.23.  The  runup  consists  of  setup  ^  ~  0.06  and  oscillatory  components  of 
amplitude  ~  0.02.  The  runup  was  not  measured  in  the  experiment  by  Cox  et  al. 
(1995). 


Table  4.6:  Primary  input  data  file  (initiated  before  breaking)  for  Cox  et  al  (1995). 
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— >  NLINES 


L13053  BDJ  June  ^95 

SBC  at  LI;  300  waves  input;  f'=.05;  and  IWAVE=3 


1 

3 

1 

1 

1 

1 

1300 

250.  300. 

490 

0.001  0.400  1.0 

.1322  2.2 

66001 

3.0  .10 

.28  0.028571 

1 


15.000 

.028571 

.05 

66001 

1 

0.5 

18 

1 

2 

3 

120 

121 

122 

180 

181 

182 

240 

241 

242 

300 

301 

302 

360 

5 

361 

362 

299.0 

299.25 

299.50 

<— ISYST 
<— IWAVE 
<— IBOT 
<— lENERG 
ITEMVA 
<— ISPAEU 
<— FINP2 
<“TSTAT,TMAX 
<— STILL 

<— DELTA ,  COURNO  ,DKAPPA 

<~-HREF,TREF 

<— NPINP 

<— APROFL,CMIXL 

<““DSEAP,SLSURF 

<— NBSEG 

<— NPOUT 
<-~NDELR 

<— NONODS 


299.75 


<— NOTIML 
300.0 


<— TIMSPA 


Table  4.7:  Primary  output  file  (initiated  before  breaking)  for  Cox  et  aL  (1995). 


L13053 

SBC  at  LI;  300  waves  input;  f^=.05;  and  IWAVE=3 


BDJ  June 


WAVE  CONDITION 


Measured  Total  Wave  Profile  at  Seaward  Boundary 


Reference  Wave  Period  = 

Reference  Wave  Height  = 

Depth  at  Seaward  Boundary  = 

Norm.  Depth  at  Seaw.  Bdr.  = 

Included  Correction  Term  CT 
0  =  no;  1  =  yes  INCLCT  =  0 

Normalized  Wave  Length  =  12.515 

"Sigma"  =  18.951 

Ursell  Number  =  73.951 

Surf  Similarity  Parameter  =  0.216 

Input  Wave  Train  from  Time=0  to  TMAX 
Computed  or  Read  at  Normalized  Rate  DELTI  = 

Parameters  of  Vertical  Velocity  Vciriations 


2.200000  sec. 
0.132200  meters 
0.280000  meters 
2.118 


0.004545 


Cubic  Profile  Parameter 

APROFL 

= 

3.000000 

Mixing  Length  Parameter 

CMIXL 

= 

0.100000 

Momentum  Flux  Coefficient 

C2 

0.548214 

Kinetic  Energy  Flux  Coeff. 

C3 

= 

-0.069420 

Energy  Dissipation  Coeff. 

CB 

= 

15.163393 

Coefficient  of  DB 

CBL 

= 

2.873676 

BOTTOM  GEOMETRY 


Norm.  Horiz.  Length  of 

Computation  Domain 
Number  of  Segments 


5.979232 

1 


SEGMENT  WBSEG(I)  TBSLOP(I)  BFFSEG(I) 

I  meters 


1  15.000000 


0.028571 


0.050000 


Table  4.7:  -  Continued 


COMPUTATION  PARAMETERS 

Normalized  DX 

0.798295D-02 

Normalized  DELTA 

= 

O.lOOOOOE-02 

Courant  Number 

= 

0.400 

Must  not  exceed  unity 

Niimerical  Damping  Coefficient 

1.0000 

Must  be  zero  or  positive 

Normalized  Computation  Duration 

TMAX  = 

300.000000 

Statistical  Calculations  Start 

when  Time  is  equal  to 

TSTAT= 

250.000000 

Total  Number  of  Spatial  Nodes 

JMAX  = 

750 

Number  of  Nodes  Along  Bottom  Below  SWL 

STILL  = 

490 

Storing  Temporal  Variations  from  Time  = 
to  TMAX  at  Normalized  Rate  DELTO  = 
Wave  Runup  Time  Series  Stored  for 

NDER  = 

Time  Series  of  ETA,  U,  and  UB 

Stored  at  NONODS  = 

Spacial  Variations  of  ETA,  U,  and  UB 


Stored  at 


NOTIML  = 


0.004545 


1  Water  Depths 


18  Nodes 


5  Time  Levels 


Maximum  time  step  = 

Minimum  time  step  = 

REFLECTION  COEFFICIENT 

ETARRMS/ETAIRMS  =  0,010 

INCIDENT  AND  REFLECTED  WAVES 


0.21941E-02 

0.13309E-02 


Inc. 

Ref. 


Max 

0.7186 

0.0098 


Min 

-0.3113 

-0.0012 


Mean 
-0 . 0277 
0.0050 


RMS 

0.3057 

0.0029 


SHORELINE  OSCILLATIONS 

Largest  Node  Number  Reached  by  Computational  Shoreline 

SMAX  =  523 


I  DELTAR(I) 
[cm] 


RUNUP (I)  RUNDOWN (I)  SETUP (I)  RMS (I) 
Ru  Rd  Zr  Rrms 


0.500 


0.099 


0.028 


0.060 


0.015 
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Figures  4.25  and  4.26  depict  the  computed  spatial  variation  of  the  normalized 
free  surface  Tj  and  the  depth  averaged  velocity  U  as  well  as  the  near  bottom  velocity 
Ub  at  five  times  throughout  the  final  wave  period.  Because  the  input  waves  specified 
at  the  seaward  boundary  are  not  exactly  periodic,  panel  one  of  Figure  4.25  exposes 
the  slight  difference  between  the  free  surfaces  at  times  t  =  299  and  t  =  300.  The 
wave  shape  changes  as  it  propagates  shoreward.  The  nonlinear  terms  in  the  shallow 
water  equations  are  responsible  for  the  steepening  of  the  wave  front  that  is  most 
apparent  in  panel  one  of  Figure  4.25.  The  magnitude  of  the  near  bottom  velocity 
shown  in  Figure  4.26  is  smaller  than  that  of  the  depth  averaged  velocity  because 
of  the  assumption  made  in  (2.44).  The  near  bottom  velocity  has  an  implausible 
discontinuity  at  the  zero  crossings  of  the  depth  averaged  velocity  caused  by  the 
abrupt  change  of  the  near  bottom  velocity  correction  based  on  (2.44a)  and  (2.44b). 
If  the  computed  momentum  flux  correction  m  does  not  approach  zero  as  the  depth 
averaged  velocity  U  approaches  zero,  a  discontinuity  will  occur  in  the  horizontal 
velocity  as  the  depth  averaged  velocity  changes  sign. 

The  comparison  of  the  computed  and  measured  time  series  of  the  normalized 
free  surface  for  the  last  ten  waves  is  depicted  in  Figure  4.27.  The  free  surface  shown 
in  panel  one  is  used  as  the  boundary  condition  for  the  computation.  Subsequent 
panels  show  the  agreement  in  phase  and  shape  at  five  measuring  lines.  Initially, 
VBREAK  overpredicts  the  phase  speed  as  seen  in  the  second  panel  at  measuring 
line  2.  After  wave  breaking,  the  model  underpredicts  the  phase  speed.  The  greatest 
accumulative  error  is  apparent  in  the  last  panel,  at  measuring  line  6.  Additionally, 
the  wave  height  is  underpredicted  at  the  measuring  lines  2-6. 

Because  the  input  wave  is  measured  data,  the  model  output  is  not  perfectly 
periodic;  comparisons  are  therefore  made  of  the  hydrodynamic  quantities  phase- 
averaged  over  the  last  50  waves.  The  comparison  of  the  phase-averaged  free  surface 
is  made  in  Figure  4.28.  Because  the  nonperiodic  components  are  relatively  small. 
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Figure  4.26:  Computed  normalized  cross-shore  depth  averaged  velocity  U  and 
near  bottom  velocity  Ub  at  5  time  levels,  t  =  299.00,  299.25,  299.50, 
299.75,  and  300:  U  (-  -  -);  Ub  (•  •  •)• 
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the  phase-averaged  free  surface  exhibits  the  same  characteristics  as  the  individual 
waves  shown  in  Figure  4.27. 

The  horizontal  velocity  data  was  taken  at  approximately  twenty  elevations 
throughout  the  water  depth.  Figures  4.29  through  4.34  compare  the  model  VBREAK 
output  with  the  phase- averaged  measured  velocity  at  three  vertical  elevations.  The 
first  panel  shows  the  comparison  near  the  wave  trough  level.  The  plot  of  the  mea¬ 
sured  phase-averaged  u  is  discontinuous  due  to  dropouts  in  the  measured  data.  Only 
the  reliable  data  near  the  trough  has  been  plotted  in  Figure  4.29  through  4.34.  The 
second  panel  of  each  figure  is  at  a  vertical  location  near  mid- depth.  Finally,  the 
last  is  a  comparison  of  the  near  bottom  velocities  at  a  distance  of  1.1  cm  above 
the  bottom.  Obviously,  the  disagreement  in  the  measured  and  computed  velocity 
variations  is  exacerbated  through  the  inaccurate  prediction  of  the  phase  speed.  For 
instance,  at  measuring  line  6,  the  velocities  are  predicted  reasonably  well  except  for 
the  mismatch  in  the  phase  between  the  computed  and  measured  data.  The  irreg¬ 
ularities  in  each  of  the  velocity  plots,  as  previously  seen  in  Sections  4.1  and  4.2, 
are  due  to  the  sign  change  in  the  bottom  velocity  correction.  The  discontinuity  in 
the  bottom  velocity  correction  accompanies  the  sign  change  in  the  depth  averaged 
velocity  U .  It  is  most  apparent  near  the  bottom  and  free  surface  where  the  values 
of  F  in  (2.33)  is  the  largest  as  shown  in  Figure  2.2. 

The  vertical  velocity  variations  at  six  phases  are  shown  in  Figures  4.35 
through  4.40.  In  each  figure,  the  first  panel  depicts  the  computed  and  measured 
horizontal  velocity  as  a  function  of  the  normalized  vertical  distance  (z  —  Zb)  above 
the  bottom.  The  computed  and  measured  free  surfaces  are  shown  as  an  empty  circle 
and  as  a  horizontal  line,  respectively.  The  measured  velocity  data  is  plotted  as  a 
solid  line,  excluding  the  data  points  that  are  deemed  unreliable  due  to  the  proximity 
of  the  free  surface.  The  second  panel  depicts  the  comparison  of  the  vertical  varia¬ 
tion  of  the  measured  vertical  velocity  with  the  value  computed  using  (2.45).  The 
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Figure  4.29:  Normalized  phase-averaged  horizontal  velocity  u  at  three  elevations 
at  measuring  line  1;  (z'  —  z^)  =  1.1,  12.1,  26.1  cm.:  Measured( — ); 
Computed  ( - ). 
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Figure  4.31:  Normalized  phase-averaged  horizontal  velocity  u  at  three  elevations 
at  measuring  line  3;  {z'  -  Zj,)  =  1.1,  8.1,  16.1  cm.:  Measured( — ); 
Computed( - ). 
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Figure  4.33:  Normalized  phase- averaged  horizontal  velocity  u  at  three  elevations 
at  measuring  line  5;  {z'  -  z[)  —  1.1,  6.1,  10.1  cm.:  Measured( — ); 
Computed  ( - ). 
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Figure  4.34:  Normalized  phase- averaged  horizontal  velocity  u  at  three  elevations 
at  measuring  line  6;  {z'  —  z^)  =  1.1,  4.1,  7.1  cm.:  Measured( — ); 


Computed (-  -  -). 


Vertical  Variations  of  Horizontal  Velocity  at  6  Phases  ml1 


Vertical  Variations  of  Vertical  Velocity  at  6  Phases  ml1 


Figure  4.35:  Vertical  variations  of  normalized  horizontal  and  vertical  velocity  at 
six  phases  at  measuring  line  1:  Measured( — );  Computed( - ). 

discrepiency  between  the  measured  and  computed  vertical  variations  are  caused,  in 
part,  by  the  aforementioned  phase  mismatch.  As  a  whole,  the  agreement  is  reason¬ 
able  in  spite  of  the  simple  vertical  velocity  profile  assumed  in  (2.33)  with  (2.40). 
This  is  probably  because  the  comparison  is  limited  below  the  wave  trough  level 
where  the  energy  dissipation  due  to  wave  breaking  is  expected  to  occur  above  the 
trough  level. 

For  completeness,  the  cross-shore  variation  of  the  maximum,  minimum  and 
mean  of  the  free  surface,  the  depth  averaged  velocity,  and  the  near  bottom  velocity 
is  depicted  in  Figure  4.41. 


Vertical  Variations  of  Horizontal  Velocity  at  6  Phases  ml2 


Vertical  Variations  of  Vertical  Velocity  at  6  Phases  ml2 


Figure  4.36:  Vertical  variations  of  normalized  horizontal  and  vertical  velocity  at 
six  phases  at  measuring  line  2:  Measured( — );  Computed( - ). 
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Vertical  Variations  of  Vertical  Velocity  at  6  Phases  ml3 


Figure  4.37:  Vertical  variations  of  normalized  horizontal  and  vertical  velocity  at 
six  phases  at  measuring  line  3:  Measured( — );  Computed( - ). 
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Vertical  Variations  of  Vertical  Velocity  at  6  Phases  ml4 


Figure  4.38:  Vertical  variations  of  normalized  horizontal  and  vertical  velocity  at 
six  phases  at  measuring  line  4:  Measured  ( — );  Computed  ( - ). 
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Figure  4.39:  Vertical  variations  of  normalized  horizontal  and  vertical  velocity  at 
six  phases  at  measuring  line  5:  Measured( — );  Computed( - ). 
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Figure  4,40:  Vertical  variations  of  normalized  horizontal  and  vertical  velocity  at 
six  phases  at  measuring  line  6:  Measured( — );  Computed( - ). 


Free  Surface  Statistics 


X 


Figure  4.41:  Computed  cross-shore  variation  of  maximum,  minimum,  and  mean 
of  the  normalized  free  surface  elevation,  depth  averaged  velocity,  and 
near  bottom  velocity. 
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Figure  4.42:  Computed  cross-shore  variation  of  normalized,  time  averaged  volume 
flux  . 

The  cross-shore  variation  of  the  time  averaged  volume  flux  is  again  used  as 
an  indicator  of  the  computational  accuracy.  Figure  4.42  demonstrates  that  this  flux 
is  zero  almost  exactly. 

The  normalized  energy  quantities  are  shown  in  Figure  4.43.  The  numerical 
dissipation  dominates  over  the  dissipation  due  to  bottom  friction  and  wave  breaking 
as  in  Section  4.2.  This  clearly  indicates  the  shortcoming  of  the  assumed  velocity 
profile  in  (2.33)  which  may  be  reasonable  below  the  trough  level  but  can  not  account 
for  the  much  larger  dissipation  occurring  above  the  trough  level. 
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Figure  4.43:  Computed  cross-shore  variation  of  normalized  specific  wave  energy, 
energy  fiux,  bottom  dissipation,  breaking  dissipation,  and  numerical 
dissipation. 


4.3.2  Model  Initiation  after  Breaking 

The  model  VBREAK  is  compared  to  the  identical  set  of  data  as  in  Sec¬ 
tion  4.3.1  except  with  the  seaward  boundary  located  at  measuring  line  3  as  depicted 
in  Figure  4.22.  The  height  of  the  incident  regular  waves  at  measuring  line  3  was 
H'  =  12.71  cm.  The  corresponding  still  water  depth  was  17.71  cm.  The  Tables  4.8 
and  4.9  represent  the  primary  input  and  output  files  from  the  model  VBREAK  . 
The  resultant  dimensionless  parameters  were:  surf  similarity  parameter  ^  =  0.22; 
the  ratio  of  horizontal  to  vertical  length  scales  a  =  19.3;  Ursell  number  =  183.0. 

Figures  4.44  through  4.60  may  be  compared  with  the  corresponding  figure 
is  Section  4.3.1.  As  a  whole,  the  agreement  between  the  measured  and  computed 
free  surface  and  velocities  is  somewhat  better  because  VBREAK  is  initiated  at  a 
measuring  line  landward  of  the  break  point  and  closer  to  measuring  lines  4,  5,  and 
6.  The  vertical  variation  of  the  horizontal  velocity  given  by  (2.33)  with  (2.40)  is 
assumed  to  be  caused  by  wave  breaking  and  it  is  hence  more  appropriate  to  take 
the  seaward  boundary  of  VBREAK  inside  the  surf  zone. 


Table  4.8:  Primary  input  data  file  (initiated  after  breaking)  for  Cox  et  ah  (1995). 
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Table  4,9:  Primary  output  file  (initiated  after  breaking)  for  Cox  et  al.  (1995). 


L33053  BDJ  June 

SBC  at  L3;  300  waves  input;  f'=,05;  and  IWAVE=3 


WAVE  CONDITION 


Measured  Total  Wave  Profile  at  Seaward  Boundary 


Reference  Wave  Period 
Reference  Wave  Height  = 

Depth  at  Seaward  Boundary  = 

Norm.  Depth  at  Seaw.  Bdr.  = 

Included  Correction  Term  CT 
0  =  no;  1  =  yes  INCLCT 

Normalized  Wave  Length 
"Sigma” 

Ursell  Number 

Surf  Similarity  Parameter 

Input  Wave  Train  from  Time=0  to 


2.200000  sec. 

0 . 127100  meters 
0.177143  meters 
1.394 


0 

15.969 

19.328 

182.968 

0.220 

TMAX 


Computed  or  Read  at  Normalized  Rate  DELTI  = 


0.004545 


Parameters  of  Vertical  Velocity  Variations 


Cubic  Profile  Parameter  APROFL  = 
Mixing  Length  Parameter  CMIXL  = 
Momentum  Flux  Coefficient  C2  = 
Kinetic  Energy  Flux  Coeff .  C3  = 
Energy  Dissipation  Coeff.  CB  = 
Coefficient  of  DB  CBL  = 


3.000000 

0.100000 

0.548214 

-0.069420 

15.163393 

2.930764 


BOTTOM  GEOMETRY 

Norm.  Horiz.  Length  of 

Computation  Domain  =  4.632540 

Number  of  Segments  =  1 


SEGMENT  WBSEG(I)  TBSLOP(I)  BFFSEG(I) 

I  meters 


1 
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0.028571 


0.050000 


Table  4.9:  -  Continued 


COMPUTATION  PARAMETERS 

Normalized  DX 

= 

0.814155D-02 

Normalized  DELTA 

= 

O.lOOOOOE-02 

C our ant  Number 

= 

0.400 

Must  not  exceed  unity 

Numerical  Damping  Coefficient 

= 

1.0000 

Must  be  zero  or  positive 

Normalized  Computation  Duration 

TMAX  = 

300.000000 

Statistical  Calculations  Start 

when  Time  is  equal  to 

TSTAT= 

250.000000 

Total  Number  of  Spatial  Nodes 

JMAX  = 

570 

Number  of  Nodes  Along  Bottom  Below  SWL 

STILL  = 

Storing  Temporal  Variations  from  Time  = 

310 

0 

to  TMAX  at  Normalized  Rate 

DELTO  = 

0.004545 

Wave  Runup  Time  Series  Stored  for 

NDER  = 

Time  Series  of  ETA,  U,  and  UB 

Stored  at  NONODS  = 

Spacial  Variations  of  ETA,  U,  and  UB 

Stored  at  NOTIML  = 


1  Water  Depths 


12  Nodes 


5  Time  Levels 


Maximum  time  step  = 

Minimum  time  step  = 

REFLECTION  COEFFICIENT 

ETARRMS/ETAIRMS  =  0.019 

INCIDENT  AND  REFLECTED  WAVES 


0.27585E-02 

0.14089E-02 


Inc. 

Ref. 


Max 
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0 . 0320 


Min 

-0 . 2849 
0.0007 


Mean 

-0.0211 

0.0169 


RMS 

0.2714 

0.0052 


SHORELINE  OSCILLATIONS 

Largest  Node  Number  Reached  by  Computational  Shoreline 

SMAX  =  364 
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Figure  4.47:  Computed  normalized  cross-shore  depth  averaged  velocity  U  and 
near  bottom  velocity  uj  at  5  time  levels,  t  =  299.00,  299.25,  299.50, 
299.7.5,  and  300:  U  (-  -  -);  Ub  (•••)• 
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Figure  4.48:  Measured  and  computed  time  series  for  the  last  ten  waves  of 
the  normalized  free  surface  at  four  measuring  lines:  Measured( — ); 


Computed!-  -  -). 
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Figure  4.51:  Normalized  phase-averaged  horizontal  velocity  u  at  three  elevations 
at  measuring  line  4;  {z'  —  z[)  =  1.1,  8.1,  13.1  cm.:  Measured( — ); 
Computed (-  -  -). 
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Vertical  Variations  of  Vertical  Velocity  at  6  Phases  ml5 


Figure  4.56:  Verti«:al  variations  of  normalized  horizontal  and  vertical  velocity  at 
six  phases  at  measuring  line  5:  Measured( — );  Computed( - ). 
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Figure  4.58:  Com])uted  cross-shore  variation  of  maximum,  minimum,  and  mean 
of  the  normalized  free  surface  elevation,  depth  averaged  velocity,  and 
near  bottom  velocity. 


Computed  cross-shore  variation  of  normalized,  time  averaged  volume 
flux  . 
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Figure  4.60:  Computed  cross-shore  variation  of  normalized  specific  wave  energy, 
energy  flux,  bottom  dissipation,  breaking  dissipation,  and  numerical 
dissipation. 


Chapter  5 


CONCLUSIONS 

The  numerical  model  VBREAK  is  developed  to  predict  the  cross-shore  and 
temporal  variations  of  the  free  surface  elevation  ?y,  the  depth-averaged  horizontal 
velocity  U,  and  the  near-bottom  horizontal  velocity  correction  ui,  associated  with  the 
momentum  flux  correction  m  due  to  the  vertical  variation  of  the  horizontal  velocity  u 
under  the  action  of  normally  incident  breaking  waves.  The  three  governing  equations 
required  for  the  computation  of  the  three  unknown  variables  are  the  depth-integrated 
continuity  and  horizontal  momentum  equations  together  with  the  new  equation  for 
the  momentum  flux  correction  m  derived  from  the  depth-integrated  wave  energy 
equation. 

The  normalized  vertical  profile  of  the  horizontal  velocity  u  outside  the  thin 
wave  boundary  layer  is  assumed  to  be  cubic  on  the  basis  of  limited  available  data. 
The  turbulent  shear  stress  outside  the  wave  boundaxy  layer  is  assumed  to  be  ex¬ 
pressed  using  the  turbulent  eddy  viscosity  whose  mixing  length  is  proportional  to 
the  instantaneous  water  depth.  Although  two  additional  empirical  parameters  are 
introduced  in  relation  to  these  assumptions,  the  computed  vertical  profiles  of  the 
horizontal  velocity  are  found  to  be  fairly  insensitive  to  these  empirical  parameters 
in  their  ranges  expected  from  limited  available  data. 

The  numerical  model  VBREAK  is  reduced  to  a  one-dimensional  model  and 
compared  to  the  previously  developed  model  IBREAK.  With  appropriate  simpli¬ 
fication  of  the  seaward  boundary  condition,  the  momentum  flux  correction  equals 


zero  identically  throughout  the  computation  domain  for  all  times.  The  computed 
results  of  the  one-dimensional  model  corresponding  to  VBREAK  are  essentially  the 
same  as  the  results  of  IBREAK  for  the  steep  slope  case.  The  comparison  of  the 
models  demonstrates  the  efficiency  and  accuracy  of  the  MacCormack  method  in  the 
solution  of  the  finite- amplitude  shallow- water  equations. 

The  model  VBREAK  is  compared  with  the  laboratory  data  presented  by 
Stive  (1980)  and  Stive  and  Wind  (1982).  The  free  surface  statistics  are  predicted 
reasonably  accurately.  The  computed  minimum  horizontal  velocities  at  several  ver¬ 
tical  locations  under  the  wave  trough  compare  well  with  the  measured  data.  The 
maximum  horizontal  velocities,  however,  are  consistently  overpredicted.  The  time- 
averaged  seaward  return  current  is  predicted  to  within  the  order  of  magnitude  only. 
Despite  the  explicit  modeling  of  energy  dissipation  due  to  wave  breaking,  most  of 
the  energy  is  dissipated  numerically.  The  assumed  velocity  profile  does  not  describe 
well  the  wave  energy  dissipation  due  to  breaking  on  a  gentle  slope.  This  is  probably 
due  to  the  absence  of  a  surface  roller  in  the  model  VBREAK  to  increase  the  energy 
dissipation  due  to  wave  breaking  and  the  shoreward  mass  flux  above  the  trough 
level. 

The  model  VBREAK  is  compared  with  the  detailed  fluid  velocity  measure¬ 
ments  of  Cox  et  al.  (1995).  VBREAK  predicts  the  general  free  surface  character¬ 
istics  well.  The  laboratory  data  of  Cox  included  phase  information  and  is  used  to 
asses  the  phase  speed  predictions  of  VBREAK.  The  computed  phase  speed  exceeds 
the  measured  speed  seaward  of  the  break  point.  Following  wave  breaking,  however, 
VBREAK  underpredicts  the  phase  speed.  The  computed  horizontal  velocities  com¬ 
pare  reasonably  well  with  the  data  except  that  differences  in  phase  generate  errors, 
especially  near  the  wave  crest.  The  agreement  with  the  data  improves  somewhat 
when  the  seaward  boundary  of  VBREAK  is  taken  to  be  landward  of  the  break 
point  because  the  vertical  horizontal  velocity  variation  in  VBREAK  is  assumed  to 


be  caused  by  wave  l^reaking.  The  discontinuity  in  the  bottom  velocity  correction 
that  accompanies  the  sign  change  in  the  depth  averaged  velocity  generates  an  unre¬ 
alistic  kink  in  the  time  series  of  the  horizontal  velocity.  As  a  whole,  the  agreement 
is  reasonable  considering  the  assumed  simple  vertical  velocity  profile. 
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